[gmx-users] Setting up an infinitely hard wall
Amit Choubey
kgp.amit at gmail.com
Thu Nov 5 02:32:08 CET 2009
Hi Everyone,
So i did set up the velocity mirror as described in the previous email.
I did few more standard calculations on water with the velocity mirror on
and compared with the past results obtained when done with other force field
& package. The results seem to be convincing. In all, the thing did work
out.
Thank you everyone and specially Berk for this.
Amit
On Fri, Oct 30, 2009 at 4:54 PM, Amit Choubey <kgp.amit at gmail.com> wrote:
> 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/20091104/c0f0a951/attachment.html>
More information about the gromacs.org_gmx-users
mailing list