I'm attempting to model a water/D2O system (298 molecules each) with varying
levels of HCl. I'm not sure exactly how to pursue this in AMBER (7.0), but
I've taken the following approach which causes sander to crash during a
minimization:
1) I've created 2 new units - H+ and Cl-. Each defines a single charged atom
(-ion species). I've attached the files generated from a SaveOff of the two
units (ions.lib)
2) A script then handles population of the (26 angstrom cubic) system
(attached - out1.pdb)
3) I've modified the parm99.dat to incorporate new values for the D and H+
atoms, and D2O bond-angle + force constants, and changed the values of the
current IM definition for the Cl- ion. I've attached (parm99.dat) for the
curious. The pertinent parts are as follows:
In the first section -
D 2.014 0.000 D in TIP3P D2O
HW 1.008 0.000 H in TIP3P water
HP 1.008 0.700 H+ ion from HCl salt
IM 35.45 2.200 assumed to be Cl- (ion minus)
In the second section defining the bond force constants:
OW-HW 553.0 0.9572 ! TIP3P water
HW-HW 553.0 1.5136 TIP3P water
OW-D 553.0 0.9572 ! TIP3P D2O
D -D 553.0 1.5136 TIP3P D2O
The third section for angles:
HW-OW-HW 100. 104.52 TIP3P water
HW-HW-OW 0. 127.74 (found in crystallographic water with 3
bonds)
D -OW-D 100. 104.52 D2O
D -D -OW 0. 127.74
And towards the end of the file:
HW 0.0000 0.0000 TIP3P water model
...
D 0.0000 0.0000 TIP3P D2O model
OW 1.7683 0.1520 TIP3P water model
4) I run leap to get a topology and the initial coordinate files. There are
only a few close-contact angle warnings when checking the unit, and a
SaveAmberParmPol does not cause anything to break and generates a topology
and coordinate file just fine.
I would appreciate any help from someone with experience in adding in ionic
species into a system as such.
The issue is that when running my script for minimization, the energy values
seem to explode and the run crashes. Below I've copied the (mdout) output
from one such run.
-------------------------------------------------------
Amber 7 SANDER Scripps/UCSF 2002
-------------------------------------------------------
| Wed Sep 13 11:41:15 2006
[-O]verwriting output
File Assignments:
| MDIN: /home/eric/MD/mdin/minimization.mdin
| MDOUT: /home/eric/MD/mdout/1.13Mhcl.min.mdout
|INPCRD: /home/eric/MD/crd/1.13Mhcl.crd
| PARM: /home/eric/MD/top/1.13Mhcl.top
|RESTRT: /home/eric/MD/crd/1.13Mhcl.min.crd
| REFC: refc
| MDVEL: mdvel
| MDEN: /home/eric/MD/en/1.13Mhcl.min.mden
| MDCRD: mdcrd
|MDINFO: mdinfo
|INPDIP: inpdip
|RSTDIP: rstdip
Here is the input file:
&cntrl
IMIN = 1,
NTX = 1,
NTPR = 200, NTWR = 200, IWRAP = 1,
NTWX = 1000, NTWV = 1000, NTWE = 1000,
NTB = 1, CUT = 8.0,
IPOL = 1,
DXM = 0.5,
MAXCYC = 80000, NCYC = 1000, NTMIN = 1, DT = 0.001
&end
--------------------------------------------------------------------------------
1. RESOURCE USE:
--------------------------------------------------------------------------------
| Flags: MPI HAS_FTN_ERFC
getting new box info from bottom of inpcrd
| peek_ewald_inpcrd: Box info found
| New format PARM file being parsed.
| Version = 1.000 Date = 09/13/06 Time = 11:32:36
NATOM = 1774 NTYPES = 5 NBONH = 1172 MBONA = 0
NTHETH = 293 MTHETA = 0 NPHIH = 0 MPHIA = 0
NHPARM = 0 NPARM = 0 NNB = 2360 NRES = 602
NBONA = 0 NTHETA = 0 NPHIA = 0 NUMBND = 2
NUMANG = 1 NPTRA = 0 NATYP = 5 NPHB = 1
IFBOX = 1 NMXRS = 3 IFCAP = 0 NEXTRA = 0
| Memory Use Allocated Used
| Real 2000000 260102
| Hollerith 400000 37858
| Integer 2000000 111503
| Max Nonbonded Pairs: 5400000
| Duplicated 0 dihedrals
| Duplicated 0 dihedrals
BOX TYPE: RECTILINEAR
--------------------------------------------------------------------------------
2. CONTROL DATA FOR THE RUN
--------------------------------------------------------------------------------
General flags:
imin = 1, nmropt = 0
Nature and format of input:
ntx = 1, irest = 0, ntrx = 1
Nature and format of output:
ntxo = 1, ntpr = 200, ntrx = 1, ntwr =
200
iwrap = 1, ntwx = 1000, ntwv = 1000,
ntwe = 1000
ioutfm = 0, ntwprt = 0, idecomp = 0,
rbornstat= 0
Potential function:
ntf = 1, ntb = 1, igb = 0,
nsnb = 25
ipol = 1, gbsa = 0
dielc = 1.00000, cut = 8.00000, intdiel = 1.00000
scnb = 2.00000, scee = 1.20000
Frozen or restrained atoms:
ibelly = 0, ntr = 0
Energy minimization:
maxcyc = 80000, ncyc = 1000, ntmin = 1
dx0 = 0.01000, dxm = 0.50000, drms = 0.00010
Polarizable options:
indmeth = 3, maxiter = 20, irstdip = 0, scaldip
= 1
diptau = 11.00000, dipmass = 0.33000
Ewald parameters:
verbose = 0, ew_type = 0, nbflag = 1, use_pme
= 1
vdwmeth = 1, eedmeth = 1, frc_int = 0,
netfrc = 1
Box X = 26.000 Box Y = 26.000 Box Z = 26.000
Alpha = 90.000 Beta = 90.000 Gamma = 90.000
NFFT1 = 30 NFFT2 = 27 NFFT3 = 27
Cutoff= 8.000 Tol = 0.100E-04
Ewald Coefficient = 0.34864
Interpolation order = 4
| PLEVEL = 1: runmd parallelization, no EKCMR
--------------------------------------------------------------------------------
3. ATOMIC COORDINATES AND VELOCITIES
--------------------------------------------------------------------------------
begin time read from input coords = 0.000 ps
Number of triangulated 3-point waters found: 0
| Atom division among processors:
| 0 1774
| Running AMBER/MPI version on 1 nodes
Sum of charges from parm topology file = -0.00000064
Forcing neutrality...
--------------------------------------------------------------------------------
4. RESULTS
--------------------------------------------------------------------------------
---------------------------------------------------
APPROXIMATING switch and d/dx switch using CUBIC SPLINE INTERPOLATION
using 5000.0 points per unit in tabled values
TESTING RELATIVE ERROR over r ranging from 0.0 to cutoff
| CHECK switch(x): max rel err = 0.2834E-14 at 2.456700
| CHECK d/dx switch(x): max rel err = 0.7554E-11 at 2.784760
---------------------------------------------------
| Local SIZE OF NONBOND LIST = 378012
| TOTAL SIZE OF NONBOND LIST = 378012
NSTEP ENERGY RMS GMAX NAME NUMBER
1 NAN NAN 0.0000E+00 O 1
BOND = 0.0559 ANGLE = 0.0002 DIHED = 0.0000
VDWAALS = nan EEL = nan HBOND =
0.0000
1-4 VDW = 0.0000 1-4 EEL = 0.0000 RESTRAINT = 0.0000
EPOLAR = nan ETHREEB = 0.0000
Dipole convergence: rms = NAN temperature = 0.00
NSTEP ENERGY RMS GMAX NAME NUMBER
200 NAN NAN 0.0000E+00 O 1
BOND = nan ANGLE = 92781.3302 DIHED = 0.0000
VDWAALS = -50.8251 EEL = nan HBOND =
0.0000
1-4 VDW = 0.0000 1-4 EEL = 0.0000 RESTRAINT = 0.0000
EPOLAR = nan ETHREEB = 0.0000
Dipole convergence: rms = 0.391E+03 temperature = 0.00
NSTEP ENERGY RMS GMAX NAME NUMBER
400 NAN NAN 0.0000E+00 O 1
BOND = nan ANGLE = 92781.3302 DIHED = 0.0000
VDWAALS = -50.8251 EEL = nan HBOND =
0.0000
1-4 VDW = 0.0000 1-4 EEL = 0.0000 RESTRAINT = 0.0000
EPOLAR = nan ETHREEB = 0.0000
Dipole convergence: rms = 0.391E+03 temperature = 0.00
NSTEP ENERGY RMS GMAX NAME NUMBER
600 NAN NAN 0.0000E+00 O 1
BOND = nan ANGLE = 92781.3302 DIHED = 0.0000
VDWAALS = -50.8251 EEL = nan HBOND =
0.0000
1-4 VDW = 0.0000 1-4 EEL = 0.0000 RESTRAINT = 0.0000
EPOLAR = nan ETHREEB = 0.0000
Dipole convergence: rms = 0.391E+03 temperature = 0.00
NSTEP ENERGY RMS GMAX NAME NUMBER
800 NAN NAN 0.0000E+00 O 1
BOND = nan ANGLE = 92781.3302 DIHED = 0.0000
VDWAALS = -50.8251 EEL = nan HBOND =
0.0000
1-4 VDW = 0.0000 1-4 EEL = 0.0000 RESTRAINT = 0.0000
EPOLAR = nan ETHREEB = 0.0000
Dipole convergence: rms = 0.391E+03 temperature = 0.00
NSTEP ENERGY RMS GMAX NAME NUMBER
1000 NAN NAN 0.0000E+00 O 1
BOND = nan ANGLE = 92781.3302 DIHED = 0.0000
VDWAALS = -50.8251 EEL = nan HBOND =
0.0000
1-4 VDW = 0.0000 1-4 EEL = 0.0000 RESTRAINT = 0.0000
EPOLAR = nan ETHREEB = 0.0000
Dipole convergence: rms = 0.391E+03 temperature = 0.00
.... RESTARTED DUE TO LINMIN FAILURE ...
.... RESTARTED DUE TO LINMIN FAILURE ...
.... RESTARTED DUE TO LINMIN FAILURE ...
.... RESTARTED DUE TO LINMIN FAILURE ...
.... RESTARTED DUE TO LINMIN FAILURE ...
FINAL RESULTS
NSTEP ENERGY RMS GMAX NAME NUMBER
1021 NAN NAN 0.0000E+00 O 1
BOND = nan ANGLE = 92781.3302 DIHED =
0.0000
VDWAALS = -50.8251 EEL = nan HBOND = 0.0000
1-4 VDW = 0.0000 1-4 EEL = 0.0000 RESTRAINT = 0.0000
EPOLAR = nan ETHREEB = 0.0000
Dipole convergence: rms = 0.391E+03 temperature = 0.00
***** REPEATED LINMIN FAILURE *****
--------------------------------------------------------------------------------
5. TIMINGS
--------------------------------------------------------------------------------
| Read coords time 0.01 ( 0.10% of Total)
| Ewald setup time 0.05 ( 3.99% of List )
| Check list validity 0.39 (33.06% of List )
| Map frac coords 0.51 (43.51% of List )
| Setup grids 0.01 ( 0.58% of List )
| Grid unit cell 0.11 ( 9.37% of List )
| Grid image cell 0.01 ( 1.27% of List )
| Build the list 0.08 ( 6.45% of List )
| Other 0.02 ( 1.77% of List )
| List time 1.17 ( 9.01% of Nonbo)
| Direct Ewald time 0.60 ( 5.07% of Ewald)
| Adjust Ewald time 0.27 ( 2.28% of Ewald)
| Self Ewald time 0.01 ( 0.09% of Ewald)
| Finish NB virial 0.12 ( 1.00% of Ewald)
| Fill Bspline coeffs 1.33 (12.55% of Recip)
| Fill charge grid 0.55 ( 5.17% of Recip)
| Scalar sum 3.02 (28.53% of Recip)
| Grad sum 0.40 ( 3.76% of Recip)
| FFT communication ti 0.01 ( 0.17% of FFT t)
| Other 5.15 (99.83% of FFT t)
| FFT time 5.16 (48.80% of Recip)
| Other 0.13 ( 1.19% of Recip)
| Recip Ewald time 10.57 (89.42% of Ewald)
| Ewald MPI wait 0.00 ( 0.03% of Ewald)
| dipole distribute ti 0.00 ( 0.03% of Ewald)
| Other 0.24 ( 2.07% of Ewald)
| Ewald time 11.82 (90.94% of Nonbo)
| Other 0.01 ( 0.06% of Nonbo)
| Nonbond force 13.00 (97.69% of Force)
| Bond energy 0.10 ( 0.72% of Force)
| Angle energy 0.13 ( 0.95% of Force)
| Dihedral energy 0.01 ( 0.06% of Force)
| FRC Collect time 0.00 ( 0.03% of Force)
| Noe calc time 0.00 ( 0.03% of Force)
| Other 0.07 ( 0.54% of Force)
| Force time 13.31 (99.29% of Runmd)
| Dipole update time 0.10 ( 0.71% of Runmd)
| Runmd Time 13.40 (96.54% of Total)
| Other 0.47 ( 3.35% of Total)
| Total time 13.88 (100.0% of ALL )
| Highest rstack allocated: 119824
| Highest istack allocated: 57876
| Setup wallclock 0 seconds
| Nonsetup wallclock 14 seconds
--
Thank you!
Eric Shamay
eshamay.uoregon.edu
Richmond Group
University of Oregon, Eugene
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
- application/octet-stream attachment: ions.li
Received on Sun Sep 17 2006 - 06:07:15 PDT
