[gmx-users] Setting up an infinitely hard wall
Amit Choubey
kgp.amit at gmail.com
Sat Oct 31 01:54:32 CET 2009
Hi Berk,
I am mailing about the issue of setting up infinite potential well.
I now want to simplify the task of inversing the velocity and mirroring the
particle by the following procedure. If a particle is found in z<delta then
i only change the z-component of the velocity v_z to mod(v_z) . This kind of
mimic's a mirror at z=0 for the velocities. This way, i wont have to worry
about the constraints because i never touch the co-ordinates of the
particles. I also think that this would work fine with pbc = xyz .
Now i present the part (do_update_md) of update.c that has to be changed.
The bold lines are the ones that i included
static void do_update_md(int start,int homenr,double dt,
t_grp_tcstat *tcstat,t_grp_acc *gstat,real nh_xi[],
rvec accel[],ivec nFreeze[],real invmass[],
unsigned short ptype[],unsigned short cFREEZE[],
unsigned short cACC[],unsigned short cTC[],
rvec x[],rvec xprime[],rvec v[],
rvec f[],matrix M,bool bExtended)
{
double imass,w_dt;
int gf=0,ga=0,gt=0;
rvec vrel;
real vn,vv,va,vb,vnrel;
real lg,xi,u;
int n,d;
if(bExtended) {
/* Update with coupling to extended ensembles, used for
* Nose-Hoover and Parrinello-Rahman coupling
* Nose-Hoover uses the reversible leap-frog integrator from
* Holian et al. Phys Rev E 52(3) : 2338, 1995
*/
for(n=start; n<start+homenr; n++) {
imass = invmass[n];
if (cFREEZE)
gf = cFREEZE[n];
if (cACC)
ga = cACC[n];
if (cTC)
gt = cTC[n];
lg = tcstat[gt].lambda;
xi = nh_xi[gt];
rvec_sub(v[n],gstat[ga].u,vrel);
for(d=0; d<DIM; d++) {
if((ptype[n] != eptVSite) && (ptype[n] != eptShell) &&
!nFreeze[gf][d]) {
vnrel = lg*(vrel[d] + dt*(imass*f[n][d] - 0.5*xi*vrel[d]
- iprod(M[d],vrel)))/(1 + 0.5*xi*dt);
/* do not scale the mean velocities u */
vn = gstat[ga].u[d] + accel[ga][d]*dt + vnrel;
v[n][d] = vn;
xprime[n][d] = x[n][d]+vn*dt;
} else {
v[n][d] = 0.0;
xprime[n][d] = x[n][d];
}
}
*if(xprime[n][2] < delta) if(v[n][2] < 0) v[n][2] = -v[n][2];*
}
} else {
/* Classic version of update, used with berendsen coupling */
for(n=start; n<start+homenr; n++) {
w_dt = invmass[n]*dt;
if (cFREEZE)
gf = cFREEZE[n];
if (cACC)
ga = cACC[n];
if (cTC)
gt = cTC[n];
lg = tcstat[gt].lambda;
for(d=0; d<DIM; d++) {
vn = v[n][d];
if((ptype[n] != eptVSite) && (ptype[n] != eptShell) &&
!nFreeze[gf][d]) {
vv = lg*(vn + f[n][d]*w_dt);
/* do not scale the mean velocities u */
u = gstat[ga].u[d];
va = vv + accel[ga][d]*dt;
vb = va + (1.0-lg)*u;
v[n][d] = vb;
xprime[n][d] = x[n][d]+vb*dt;
} else {
v[n][d] = 0.0;
xprime[n][d] = x[n][d];
}
}
*if(xprime[n][2] < delta) if(v[n][2] < 0) v[n][2] = -v[n][2];*
}
}
}
I assume that xprime[n][2] and v[n][2] refer to the z-components of position
and velocity. xprime means the new co-ordinates. Do you think this change is
enough to implement the stuff i was talking about? I must acknowledge that i
have very little knowledge about the gromacs codes and my views might be
very short sighted. So, please do let me know any of the shortcomings that
you notice.
If everything is fine with the above, could you also suggest me how to
recompile everything.
I cannot thank you all enough for helping me out.
Amit
On Fri, Oct 23, 2009 at 1:09 AM, Berk Hess <gmx3 at hotmail.com> wrote:
> Hi,
>
> If you have a wall potential, you can reweight the configurations
> with weight 0 if one or more particles are beyond the wall
> and exp(Vwall/kT) if no particles are beyond the wall.
> This will only work efficiently if you choose the wall potential in a smart
> way.
> I managed to get an efficiency of 70%.
>
> If you have constraints, an atom that goes beyond the wall is directly
> coupled to the atoms it is constrained to. You will have to work out
> the equations for such a situation. It might be as simple that you can
> just correct x and v for the atom that went beyond the wall and then
> apply the constraints as normal.
>
> Berk
>
> ------------------------------
> Date: Fri, 23 Oct 2009 00:26:08 -0700
> Subject: Re: [gmx-users] Setting up an infinitely hard wall
>
> From: kgp.amit at gmail.com
> To: gmx-users at gromacs.org
>
> Hi Berk,
>
> Thank you for the response.
>
> I have obtained distributions with in infinite wall by simulating with an
> softer wall
> and unbiasing with configurations with the wall potential.
>
>
> Could you explain the above ? I am not sure if i get your point.
>
>
> But if you need the dynamics, or you want less hassle during the
> simulation,
> you will have to do a bit more effort in coding an infinitely hard wall.
>
> You will not only have to inverse the velocity, but also mirror the
> position
> of the particle in the wall.
>
>
> That's true. I am sorry i didnt mention that.
>
>
> This should on require a few lines of code in do_update_md
> in src/mdlib/update.c.
>
>
> Great, I will try this as soon as possible.
>
>
> Note that this will only work easily when you do not have constraints
> present.
> With constraints things get much more complicated.
>
>
> I do have constraints present. Thank you for pointing this. I am working
> with SPC water and I should, in principle be able to figure out the
> co-ordinates of all atoms in the molecule if I am given one of the water's
> atom's co-ordinate. Does that sound ok?
>
> Thanks for the input again,
> Amit
>
>
>
> Berk
>
> ------------------------------
> Date: Thu, 22 Oct 2009 14:34:06 -0700
> From: kgp.amit at gmail.com
> To: gmx-users at gromacs.org
> Subject: [gmx-users] Setting up an infinitely hard wall
>
>
> Hi everyone,
>
> I am sending this email again hoping for any quick input for my question.
>
>
>
> I have been trying to set up an "infinitely" hard potential wall.
>
> I tried to use the available wall options and could not really get it to do
> what i needed. I wanted a steep repulsive potential but when i created that,
> the system was blowing up, reason being that it requires smaller time step
> and i cant afford to have smaller time step.
>
> My idea to overcome this issue is to just reverse the velocity of the
> particle whenever it hits the wall. I am not sure if there is any thing in
> GROMACS which does this but any suggestions will be very helpful.
>
> If there is nothing set up for something like above, i would like to play
> around with the code to figure it out. If this is the case, could somebody
> direct me to a starting point.
>
> Thank you
> Amit
>
>
> ------------------------------
> New Windows 7: Find the right PC for you. Learn more.<http://windows.microsoft.com/shop>
>
> _______________________________________________
> gmx-users mailing list gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>
>
>
> ------------------------------
> Express yourself instantly with MSN Messenger! MSN Messenger<http://clk.atdmt.com/AVE/go/onm00200471ave/direct/01/>
>
> _______________________________________________
> gmx-users mailing list gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20091030/8840ccc2/attachment.html>
More information about the gromacs.org_gmx-users
mailing list