# [gmx-users] Re: Re: Re: Definition of "polarization"

Yinghong xieyh at hkusua.hku.hk
Fri Dec 16 21:37:17 CET 2005

```David:

>
> /
> />/  >/ Dear Gmx'ers:
> />/ />/
> />/ />/ As mentioned in the publication about shell water model, a shell
> />/ />/ particle is connected to a dummy atom by a spring-like
> connection with
> />/ />/ the following relationship:
> />/ />/
> />/ />/ 1. K=sqr(qS)/(4*PHI*Epsilon*alpha),
> />/ />/ 2. rsd=(4*PHI*Epsilon*alpha)*E / qS;  (rsd is the distance between
> />/ shell
> />/ />/ and dummy particles at any moment, and E is refered to
> electrical field/
>
>  >/ />/ strength).
> />/ />/
> />/ />/ Anyone can tell me where is the definition of two equations above,
> />/ />/ especially in the case of* "isotropic polarization",* in GMX source
> />/ code.
> />/ />/ Because in my simulation, what I want to polarized is not a water
> />/ />/ molecule. So I need a definite instruction, any suggestion?
> />/ /
> />/  >gmx/src/gmxlib/bondfree.c
> />/
> />/ In this file, I can only locate the definition of K, and where is
> />/ the formula for calculating rsd? Pls see the following code, it is
> />/ refered to the case of isotropic polarization, ok?
> /rsd is the input for the equation, as it is just the instantaneous
> distance between shell and particle.
>  >/
> />/
> />/ />/ Secondly, how can I read the calculated value of spring constant
> K by
> />/ />/ GMX after I defined "alpha" and "qS"? Can "debug" realize it? If
> yes,
> />/ />/ how? Although I have used GMX for a long time, but that is the first
> />/ />/ time for me to encounter such a problem.
> />/ /
> />/  >gmxdump -s topol.tpr | less
> />/  >search for POL
> />/ Yes, I can find the defined value for "polarization: alpha", is it
> />/ possible for me to know the calculated value of spring constant K? If
> />/ yes, how?
> /
>
>  > It is not printed anywhere. If you want to see it you have to turn on
>  > the debug flag for mdrun.
>
> Very sorry to bother you again, could you point it out much more
> detailed? I have never used the "-debug" flag.
>
> For example, in original, I used "mdrun -v -s full -e full -o full -c
> after_full -g full" to run my simulation.
>
> Where should "-debug" flag put in? And, where or which file can I find
> the calculated value of spring constant K?

> mdrun -debug (other options)
> you will get a file mdrun.log contaning lots of debug information,
> search for ksh.

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:

real polarize(int nbonds,
t_iatom forceatoms[],t_iparams forceparams[],
rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
matrix box,real lambda,real *dvdlambda,
t_mdatoms *md,int ngrp,real egnb[],real egcoul[],
t_fcdata *fcd)
{
int  i,m,ki,ai,aj,type;
real dr,dr2,fbond,vbond,fij,vtot,ksh;
rvec dx;
ivec dt;

vtot = 0.0;
for(i=0; (i<nbonds); ) {
type = forceatoms[i++];
ai   = forceatoms[i++];
aj   = forceatoms[i++];
ksh  = sqr(md->chargeT[aj])*ONE_4PI_EPS0/forceparams[type].polarize.alpha;

ki   = pbc_rvec_sub(x[ai],x[aj],dx); /*   3   */
dr2  = iprod(dx,dx);   /*   5  */
dr   = dr2*invsqrt(dr2);          /*  10  */

*dvdlambda += harmonic(ksh,ksh,0,0,dr,lambda,&vbond,&fbond);  /*  19  */

if (dr2 == 0.0)
continue;

vtot  += vbond;/* 1*/
fbond *= invsqrt(dr2);   /*   6  */

if (g) {
ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
ki=IVEC2IS(dt);
}
for (m=0; (m<DIM); m++) {   /*  15  */
fij=fbond*dx[m];
f[ai][m]+=fij;
f[aj][m]-=fij;
fr->fshift[ki][m]+=fij;
fr->fshift[CENTRAL][m]-=fij;
}
}     /* 59 TOTAL */
return vtot;
}

Thanks and regards,
Xie Yinghong
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20051216/11d1d214/attachment.html>
```