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

[David van der Spoel]

Yinghong wrote:
> 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;

Weird, my 3.3 code has this extra line, please
add this line here and recompile:

     if (debug)
       fprintf(debug,"POL: ai + %d aj = %d ksh = %.3f\n",ai,aj,ksh);


>  
>     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
> 
>  
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> gmx-users mailing list
> gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-users-request at gromacs.org.


-- 
David.
________________________________________________________________________
David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group,
Dept. of Cell and Molecular Biology, Uppsala University.
Husargatan 3, Box 596,  	75124 Uppsala, Sweden
phone:	46 18 471 4205		fax: 46 18 511 755
spoel at xray.bmc.uu.se	spoel at gromacs.org   http://xray.bmc.uu.se/~spoel
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++