<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML><HEAD>
<META http-equiv=Content-Type content="text/html; charset=gb2312">
<META content="MSHTML 6.00.2900.2802" name=GENERATOR></HEAD>
<BODY>
<DIV>David:</DIV>
<DIV><BR>><I> <BR></I>><I> / <BR></I>><I> />/ >/ Dear
Gmx'ers:<BR></I>><I> />/ />/<BR></I>><I> />/ />/ As mentioned
in the publication about shell water model, a shell<BR></I>><I> />/ />/
particle is connected to a dummy atom by a spring-like <BR></I>><I>
connection with<BR></I>><I> />/ />/ the following
relationship:<BR></I>><I> />/ />/<BR></I>><I> />/ />/ 1.
K=sqr(qS)/(4*PHI*Epsilon*alpha),<BR></I>><I> />/ />/ 2.
rsd=(4*PHI*Epsilon*alpha)*E / qS; (rsd is the distance
between<BR></I>><I> />/ shell<BR></I>><I> />/ />/ and dummy
particles at any moment, and E is refered to <BR></I>><I> electrical
field/<BR></I>><I> <BR></I>><I> >/ />/
strength).<BR></I>><I> />/ />/<BR></I>><I> />/ />/ Anyone can
tell me where is the definition of two equations above,<BR></I>><I> />/
/>/ especially in the case of* "isotropic polarization",* in GMX source</DIV>
<DIV></I>><I> />/ code.<BR></I>><I> />/ />/ Because in my
simulation, what I want to polarized is not a water<BR></I>><I> />/ />/
molecule. So I need a definite instruction, any suggestion?<BR></I>><I>
/>/ /<BR></I>><I> />/
>gmx/src/gmxlib/bondfree.c<BR></I>><I> />/ <BR></I>><I> />/ In
this file, I can only locate the definition of K, and where is<BR></I>><I>
/>/ the formula for calculating rsd? Pls see the following code, it
is<BR></I>><I> />/ refered to the case of isotropic polarization,
ok?<BR></I>><I> /rsd is the input for the equation, as it is just the
instantaneous<BR></I>><I> distance between shell and
particle.<BR></I>><I> >/ <BR></I>><I> />/ <BR></I>><I>
/>/ />/ Secondly, how can I read the calculated value of spring constant
<BR></I>><I> K by<BR></I>><I> />/ />/ GMX after I defined "alpha"
and "qS"? Can "debug" realize it? If <BR></I>><I> yes,<BR></I>><I> />/
/>/ how? Although I have used GMX for a long time, but that is the first
</I>><I> />/ />/ time for me to encounter such a
problem.<BR></I>><I> />/ /<BR></I>><I> />/ >gmxdump -s
topol.tpr | less<BR></I>><I> />/ >search for POL<BR></I>><I>
/>/ Yes, I can find the defined value for "polarization: alpha", is
it<BR></I>><I> />/ possible for me to know the calculated value of spring
constant K? If<BR></I>><I> />/ yes, how?<BR></I>><I> /<BR></I>><I>
<BR></I>><I> > It is not printed anywhere. If you want to see it you
have to turn on<BR></I>><I> > the debug flag for
mdrun.<BR></I>><I> <BR></I>><I> Very sorry to bother you again, could you
point it out much more <BR></I>><I> detailed? I have never used the "-debug"
flag.<BR></I>><I> <BR></I>><I> For example, in original, I used "mdrun -v
-s full -e full -o full -c <BR></I>><I> after_full -g full" to run my
simulation.<BR></I>><I> <BR></I>><I> Where should "-debug" flag put in?
And, where or which file can I find <BR></I>><I> the calculated value of
spring constant K?<BR></I></DIV>
<DIV>> mdrun -debug (other options)<BR>> you will get a file mdrun.log
contaning lots of debug information, <BR>> search for ksh.</DIV>
<DIV> </DIV>
<DIV>I have tried "mdrun -debug", but I can not find out "ksh" item in
mdrun.log. Because what I adopt is "isotropic polarization", maybe, there is no
output for "ksh" in mdrun.log? At least, I can not find relevant
information about the "ksh" output while debugging in the
following code:</DIV>
<DIV> </DIV>
<DIV> real polarize(int nbonds,<BR>
t_iatom forceatoms[],t_iparams
forceparams[],<BR> rvec x[],rvec
f[],t_forcerec *fr,t_graph *g,<BR> matrix
box,real lambda,real *dvdlambda,<BR>
t_mdatoms *md,int ngrp,real egnb[],real
egcoul[],<BR> t_fcdata *fcd)<BR>{<BR>
int i,m,ki,ai,aj,type;<BR> real
dr,dr2,fbond,vbond,fij,vtot,ksh;<BR> rvec dx;<BR> ivec dt;</DIV>
<DIV> </DIV>
<DIV> vtot = 0.0;<BR> for(i=0; (i<nbonds); )
{<BR> type = forceatoms[i++];<BR>
ai = forceatoms[i++];<BR> aj =
forceatoms[i++];<BR> ksh =
sqr(md->chargeT[aj])*ONE_4PI_EPS0/forceparams[type].polarize.alpha;<BR>
<BR> ki =
pbc_rvec_sub(x[ai],x[aj],dx); /* 3
*/<BR> dr2 =
iprod(dx,dx); /*
5 */<BR> dr =
dr2*invsqrt(dr2); /*
10 */</DIV>
<DIV> </DIV>
<DIV> *dvdlambda +=
harmonic(ksh,ksh,0,0,dr,lambda,&vbond,&fbond); /* 19
*/</DIV>
<DIV> </DIV>
<DIV> if (dr2 == 0.0)<BR>
continue;<BR> <BR> vtot += vbond;/*
1*/<BR> fbond *= invsqrt(dr2); /*
6 */</DIV>
<DIV> </DIV>
<DIV> if (g) {<BR>
ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);<BR>
ki=IVEC2IS(dt);<BR> }<BR> for (m=0;
(m<DIM); m++) { /*
15 */<BR>
fij=fbond*dx[m];<BR>
f[ai][m]+=fij;<BR>
f[aj][m]-=fij;<BR>
fr->fshift[ki][m]+=fij;<BR>
fr->fshift[CENTRAL][m]-=fij;<BR> }<BR>
} /* 59 TOTAL */<BR> return
vtot;<BR>}</DIV>
<DIV> </DIV>
<DIV> </DIV>
<DIV> </DIV>
<DIV>Thanks and regards,</DIV>
<DIV>Xie Yinghong</DIV>
<DIV><BR> </DIV></BODY></HTML>