|
Thank for your help, all the suggestions and explanations . I will have to check all of them. I’ll let you know about the results. all best I am not sure what's going on;
it could be any number of things.
However, I would probably start by trying to create a topology which will work for free_energy=no (which corresponds to lambda=0 which corresponds to no perturbation (at least in the case when soft core is not used)). If your topology works properly then (which it should), then there are a couple of other things I would suggest. 1) sc-alpha IMO is usually better at 0.5 than 1.5 for particle annhilation, at least with sc-power=1.0 (the default). I've done some checking on this and this seems to generally be the optimal choice of parameters (which is also what Michael Shirts has found). You might check if changing this helps. (In my experience sc-alpha=1.5 can lead to really steep slopes in dv/dlambda near lambda=0 and lambda=1, which could lead to large forces and some of the instabilities you seem to be seeing). 2) There are a number of people who do exactly what you're doing and turn off charges and van der Waals interactions at the same time, but I don't think it's really correct to do so. First, you want a smooth dv/dlambda curve, but if you adjust soft core parameters to give you a relatively smooth curve for turning off van der Waals interactions, it gives you a very un-smooth curve for turning off charge interactions, meaning that if you do both of them at the same time, your curve won't be very smooth. In principle, there can be worse problems, too: Soft core allows atoms to start interpenetrating somewhere before lambda=1, when the charges are still on to some extent. So now you have a Cl- with very weak vdW interactions such that other atoms can overlap with it -- atoms like, for example, the hydrogens on water (which have the opposite charge and typically no vdW radii either). So now you're putting a positively charged atom almost on top of a (somewhat) negatively charged ion, and bad things happen (usually the simulation blows up. You could describe this as sort of a fusion event... :) ). My typical approach is to turn off charges first in a separate calculation without using soft core. The electrostatic interactions are usually a very smooth function of lambda, so this can usually be done with only 5 lambda values or so. THEN use soft core to turn off the vdW interactions. Again, there are two advantages to this: (a) Although it IS harder to run two sets of calculations instead of one, you can usually get by with fewer lambda values on each since now both transformations are relatively smooth functions of lambda. So it's two sets of calculations but you may be able to do them with the same total number of lambda values. (b) You avoid charge-charge overlap. IMO your problem is more likely to be the first than the second, since you are seeing problems at lambda=0 as well. Do let me know what you find out, though, especially if it is the second problem (I need to demonstrate the charge-charge overlap problem, but it is only sporadic in everything I've tried, so I'd like to find a system where it's reproducible). Incidentally, are you also using PME for long range electrostatics? Best wishes, David Mobley UCSF On 3/30/06, P <pyt1 at op.pl> wrote: > > Hi all. > I'm trying FEP In gromacs3.3. My system (3nm; 3nm; 3nm) consists > of one Na, Cl, DUM ion and solvent molecules. > During my free energy calculations I want to perturb: > 1) Cl- into a DUM atom > 2) DUM into Cl- > > > I'm using ffgmx. I've added new atom DUM into ffgmx.atp: > ….. > SD 32.06000 ; DMSO Sulphur > OD 15.99940 ; DMSO Oxygen > CD 15.03500 ; DMSO Carbon > CHE 12.01100 ; HEME RING CARBON > MNH3 0 ; Dummy mass in rigid tetraedrical NH3 group > MW 0 ; Dummy mass in rigid tyrosine rings > IW 0 ; Dummy particle in TIP4P etc. > DUM 0 ; moj atom > ….. > > I've updated ffgmxnb.itp: > …. > OWT3 15.99940 0.000 A 0.24889E-02 0.24352E-05 > SD 32.06000 0.000 A 0.10561E-01 0.21499E-04 > OD 15.99940 0.000 A 0.22715E-02 0.75147E-06 > CD 15.03500 0.000 A 0.90507E-02 0.21758E-04 > CHE 12.01100 0.000 A 0.23402E-02 0.33740E-05 > MNH3 0 0.000 A 0.0 0.0 > MW 0 0.000 A 0.0 0.0 > DUM 0 0.000 A 0.0 0.0 > …. > > My dual topology for Cl and DUM atoms: > > [ moleculetype ] > ; molname nrexcl > Cl 1 > > [ atoms ] > ; id at type res nr residu name at name cg nr charge > 1 DUM 1 DUM DUM 1 0 35.45300 > Cl -1 35.45300 > > > … and > > [ moleculetype ] > ; molname nrexcl > ClB 1 > > [ atoms ] > ; id at type res nr residu name at name cg nr charge > 1 Cl 1 Cl Cl 1 -1 35.45300 > DUM 0 35.45300 > > > I did calculations for 20 lambda values 0.05, 0.10 , 0.15 …. > 0.95 and for those values delta free energy looks reasonable. > My problem is that I get error every time I try to run > FEP for lambda = 1 or lambnda = 0. > Am I doing something wrong. Please give me some advice. > (I'm using PME, DUM and Cl atom are restrained, step 1fs) > > > My ERROR message for lambda = 0 is : > … > …. > > Large VCM(group Cl): -38.08457, 8.92837, -19.97197, ekin-cm: > 6.83899e+04 > Large VCM(group Na): 0.04709, -0.21100, -0.29647, ekin-cm: > 1.54762e+00 > > t = 0.002 ps: Water molecule starting at atom 61 can not be settled. > Check for bad contacts and/or reduce the timestep.Wrote pdb files with > previous and current coordinates > Large VCM(group SOL): -0.01922, -0.00999, -0.02540, ekin-cm: > 8.98334e+00 > Large VCM(group Cl): 4.35556, 2.33494, 5.86521, ekin-cm: > 2.08547e+03 > Large VCM(group Na): 0.04590, -0.19717, -0.27716, ekin-cm: > 1.35412e+00 > Large VCM(group SOL): -0.05937, 0.00317, -0.07617, ekin-cm: > 7.52736e+01 > Large VCM(group Cl): 13.48634, -0.66233, 17.40376, ekin-cm: > 1.72022e+04 > Large VCM(group Na): 0.04470, -0.18096, -0.25439, ekin-cm: > 1.14329e+00 > Large VCM(group SOL): -0.00393, -0.03332, 0.01434, ekin-cm: > 1.07341e+01 > Large VCM(group Cl): 0.88076, 7.63229, -3.18263, ekin-cm: > 2.45181e+03 > > t = 0.003 ps: Water molecule starting at atom 631 can not be settled. > Check for bad contacts and/or reduce the timestep.Wrote pdb > files with previous and current coordinates > Large VCM(group Cl): 6.73354, 14.12955, -14.62614, ekin-cm: > 1.62697e+04 > > > Thank you for your help > ALL BEST > |
_______________________________________________ gmx-users mailing list gmx-users@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 [EMAIL PROTECTED] Can't post? Read http://www.gromacs.org/mailing_lists/users.php
