Hi, For the dynamics to work you indeed need smaller forces (on the order of 10^3 in GMX units). Using flexible water molecules I was able to get this for your NaCl model. Just add: define = -DFLEXIBLE To your input.
This should work also to fill in the voids I guess. Ran -- ------------------------------------------------------ Ran Friedman Postdoctoral Fellow Computational Structural Biology Group (A. Caflisch) Department of Biochemistry University of Zurich Winterthurerstrasse 190 CH-8057 Zurich, Switzerland Tel. +41-44-6355559 Email: [email protected] Skype: ran.friedman ------------------------------------------------------ Justin A. Lemkul wrote: > > Use a different box size. I replicated your problem, but your run > completes successfully with a box set up with editconf -d 1 instead of > -d 3.3. > > If you set nstxout = 1 during the EM process, you'll see the > problematic water molecule become unstable. It looks as if there is a > small void in the solvent (due to the way genbox tries to stack > solvent configurations) and your water molecule simply can't find a > good orientation within that void. > > -Justin > > Ehud Schreiber wrote: >> Dear GROMACS users, >> >> >> >> When trying to simulate a pair of interacting proteins in water, I >> have encountered problems that ultimately resulted in the simulation >> crashing. I then tried to simplify the system as far as possible >> while retaining the problem; I now believe the problem (or at least a >> part of it) lies in the energy minimization step (the first molecular >> dynamics one). Specifically, the forces encountered during this step >> are very large, and some water molecules (which are supposed to be >> rigid) become giant and misshapen. >> >> >> >> In more details: >> >> 1) I use GROMACS version 4.0.7, single precision, on a server with >> two Intel x86-64 processors and the redhat 5.4 linux OS. >> >> 2) I created a PDB file, called NaCl.pdb, with only two “atoms”, >> actually Na+ and Cl- ions separated by their distance in the lattice >> of salt: >> >> >> >> HET NA A 1 1 >> >> HET CL A 1 1 >> >> HETNAM NA SODIUM ION >> >> HETNAM CL CHLORIDE ION >> >> FORMUL 1 NA NA 1+ >> >> FORMUL 2 CL CL 1- >> >> HETATM 1 NA NA A 1 -1.41 0.0 0.0 1.00 >> 0.0 NA >> >> HETATM 2 CL CL A 1 +1.41 0.0 0.0 1.00 >> 0.0 CL >> >> >> >> 3) I use the tip3p water model: >> >> >> >> pdb2gmx -f NaCl.pdb -water tip3p >> >> >> >> 4) I create the box: >> >> >> >> editconf -f conf.gro -bt dodecahedron -d 3.3 -o box.gro >> >> >> >> 5) I add water using spc216, creating the saltwater.gro file (which >> seems O.K. by inspection): >> >> >> >> genbox -cp box.gro -cs spc216.gro -p topol.top -o saltwater.gro >> >> >> >> 6) I create the energy minimization parameter file em.mdp: >> >> >> >> ------em.mdp------ >> >> integrator = steep >> >> nsteps = 200 >> >> nstlist = 10 >> >> rlist = 1.0 >> >> coulombtype = pme >> >> rcoulomb = 1.0 >> >> vdwtype = Cut-off >> >> rvdw = 1.0 >> >> nstenergy = 10 >> >> ------------------ >> >> >> >> 7) I prepare the em.tpr file for the energy minimization run: >> >> >> >> grompp -f em.mdp -p topol.top -c saltwater.gro -o em.tpr >> >> >> >> 8) I run the energy minimization step: >> >> >> >> mdrun -v -deffnm em >> >> >> >> 9) Looking at the em.log file I see that this step converged to >> machine precision but did not have maximal force < 10: >> >> >> >> … >> >> Enabling SPC water optimization for 7561 molecules. >> >> … >> >> Max number of connections per atom is 2 >> >> Total number of connections is 30244 >> >> Max number of graph edges per atom is 2 >> >> Total number of graph edges is 30244 >> >> Going to use C-settle (7561 waters) >> >> wo = 0.333333, wh =0.333333, wohh = 3, rc = 0.075695, ra = 0.0390588 >> >> rb = 0.0195294, rc2 = 0.15139, rone = 1, dHH = 0.15139, dOH = 0.09572 >> >> … >> >> Stepsize too small, or no change in energy. >> >> Converged to machine precision, >> >> but not to the requested precision Fmax < 10 >> >> … >> >> Steepest Descents converged to machine precision in 36 steps, >> >> but did not reach the requested Fmax < 10. >> >> Potential Energy = -3.4678925e+05 >> >> Maximum force = 6.4623531e+05 on atom 11052 >> >> Norm of force = 5.4643726e+03 >> >> >> >> 10) Looking at the em.gro file I see one monstrous water molecule >> (no. 3686); e.g., it has |HW2-OW| = 3.384876 nm, while the normal >> distance is about 0.1 nm. Its HW2 atom (no. 11054) is close to >> another water molecule (no. 5849), e.g., 0.047 nm from the latter’s >> HW2 atom (no. 17543): >> >> >> >> … >> >> 3686SOL OW11052 4.348 3.778 -0.629 >> >> 3686SOL HW111053 5.360 2.601 0.505 >> >> 3686SOL HW211054 6.518 1.650 0.861 >> … >> >> 5849SOL OW17541 6.525 1.698 0.900 >> >> 5849SOL HW117542 6.606 1.649 0.918 >> >> 5849SOL HW217543 6.481 1.648 0.832 >> >> … >> >> >> >> 11) During the simulation, several files called stepnnl.pdb were >> produced for problematic steps, where nn=11,15,19 and l=b,c. For >> example, the file step19c.pdb indeed shows a problematic water >> molecule no. 3686, while step19b.pdb does not. Likewise, the earlier >> step11c.pdb shows a problematic water molecule no. 3266 while >> step11b.pdb seems proper. The stdout/stderr of mdrun contains >> warnings like the following: >> >> >> >> … >> >> t = 0.019 ps: Water molecule starting at atom 11052 can not be settled. >> >> Check for bad contacts and/or reduce the timestep. >> >> … >> >> >> >> 12) Reducing the neighbor list update time, i.e., setting nstlist = >> 1, does not produce any change. >> >> >> >> 13) Trying to use the conjugate gradient integrator instead of >> steepest descent, i.e., setting integrator = cg, is even worse - the >> running crashes: >> >> >> >> … >> >> Step 16, Epot=-1.258771e+35, Fnorm= nan, Fmax= inf (atom >> 14493) >> >> Segmentation fault >> >> Exit 139 >> >> >> >> >> >> >> >> So, am I doing something wrong? How can I avoid these problems? >> >> >> >> Thanks, >> >> Ehud Schreiber. -- gmx-users mailing list [email protected] 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 [email protected]. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
