if you look at the outputs, the temperature rises before you get the NaN
errors. Why the temperature isn't stable is something you'll need to check.
Is the PDB file from experiments, with no modified force field parameters?
an unreliable pdb or non-standard parameters can cause issues. you also
might want to look at the potential energy components when the temperature
starts to rise- have any of the energies increased?
On Thu, Feb 9, 2023 at 10:26 AM Prithviraj Nandigrami <
prithviraj.nandigrami.gmail.com> wrote:
> I was able to resolve the issue of LINMIN failure at that minimization
> stage and move on to the subsequent steps of relaxation protocol.
>
> I faced another error at the relaxation stage when reducing the restraint
> on the backbone. The input file 6md.in looks like:
>
> &cntrl
> imin = 0, nstlim = 1000000, dt = 0.001,
> irest = 0, ntx = 1, ig = -1,
> tempi = 300.0, temp0 = 300.0,
> ntc = 2, ntf = 2, tol = 0.00001,
> ntwx = 10000, ntwe = 0, ntwr = 1000, ntpr = 10,
> cut = 8.0, iwrap = 0,
> ntt =3, gamma_ln=1.0, ntb = 2, ntp = 1,
> nscm = 0, barostat = 2,
> ntr=1,
> restraintmask=".OP1,OP2,OP3,P,C1',C2',C3',C4',C5',O2',O3',O4',O5'",
> *restraint_wt=10.0*,
> ioutfm=1, ntxo=2,
> /
>
> This stage works fine. However, in the next stage I reduce the value of
> restraint_wt on the backbone as follows:
>
> &cntrl
> imin = 0, nstlim = 1000000, dt = 0.001,
> irest = 1, ntx = 5, ig = -1,
> temp0 = 300.0,
> ntc = 2, ntf = 2, tol = 0.00001,
> ntwx = 10000, ntwe = 0, ntwr = 1000, ntpr = 10,
> cut = 8.0, iwrap = 0,
> ntt =3, gamma_ln=1.0, ntb = 2, ntp = 1,
> nscm = 0, barostat = 2,
> ntr=1,
> restraintmask=".OP1,OP2,OP3,P,C1',C2',C3',C4',C5',O2',O3',O4',O5'",
> *restraint_wt=1.0*,
> ioutfm=1, ntxo=2,
> /
>
> This step runs for a while and I checked the Density and vlimit values,
> they look fine. However, during the 7md stage, the temperature values lokk
> like the following for the last few steps after which the simulation
> crashes:
>
>
> 7md.out: NSTEP = 61090 TIME(PS) = 1061.090 TEMP(K) = 300.65
> PRESS = 0.0
> 7md.out: NSTEP = 61100 TIME(PS) = 1061.100 TEMP(K) = 298.95
> PRESS = 0.0
> 7md.out: NSTEP = 61110 TIME(PS) = 1061.110 TEMP(K) = 299.51
> PRESS = 0.0
> 7md.out: NSTEP = 61120 TIME(PS) = 1061.120 TEMP(K) = 300.00
> PRESS = 0.0
> 7md.out: NSTEP = 61130 TIME(PS) = 1061.130 TEMP(K) = 300.42
> PRESS = 0.0
> 7md.out: NSTEP = 61140 TIME(PS) = 1061.140 TEMP(K) = 358.81
> PRESS = 0.0
> 7md.out: NSTEP = 61150 TIME(PS) = 1061.150 TEMP(K) = 376.60
> PRESS = 0.0
> 7md.out: NSTEP = 61160 TIME(PS) = 1061.160 TEMP(K) = 400.54
> PRESS = 0.0
> 7md.out: NSTEP = 61170 TIME(PS) = 1061.170 TEMP(K) = 423.08
> PRESS = 0.0
> 7md.out: NSTEP = 61180 TIME(PS) = 1061.180 TEMP(K) = 417.14
> PRESS = 0.0
> 7md.out: NSTEP = 61190 TIME(PS) = 1061.190 TEMP(K) = 423.73
> PRESS = 0.0
> 7md.out: NSTEP = 61200 TIME(PS) = 1061.200 TEMP(K) = 407.54
> PRESS = 0.0
> 7md.out: NSTEP = 61210 TIME(PS) = 1061.210 TEMP(K) = 388.32
> PRESS = 0.0
> 7md.out: NSTEP = 61220 TIME(PS) = 1061.220 TEMP(K) = 365.34
> PRESS = 0.0
> 7md.out: NSTEP = 61230 TIME(PS) = 1061.230 TEMP(K) = 340.88
> PRESS = 0.0
> 7md.out: NSTEP = 61240 TIME(PS) = 1061.240 TEMP(K) = 324.64
> PRESS = 0.0
> 7md.out: NSTEP = 61250 TIME(PS) = 1061.250 TEMP(K) = 317.69
> PRESS = 0.0
> 7md.out: NSTEP = 61260 TIME(PS) = 1061.260 TEMP(K) = 321.64
> PRESS = 0.0
> 7md.out: NSTEP = 61270 TIME(PS) = 1061.270 TEMP(K) = 327.02
> PRESS = 0.0
> 7md.out: NSTEP = 61280 TIME(PS) = 1061.280 TEMP(K) = 339.70
> PRESS = 0.0
> 7md.out: NSTEP = 61290 TIME(PS) = 1061.290 TEMP(K) = 354.73
> PRESS = 0.0
> 7md.out: NSTEP = 61300 TIME(PS) = 1061.300 TEMP(K) = 366.19
> PRESS = 0.0
> 7md.out: NSTEP = 61310 TIME(PS) = 1061.310 TEMP(K) = 376.26
> PRESS = 0.0
> 7md.out: NSTEP = 61320 TIME(PS) = 1061.320 TEMP(K) = 386.98
> PRESS = 0.0
> 7md.out: NSTEP = 61330 TIME(PS) = 1061.330 TEMP(K) = 390.69
> PRESS = 0.0
> 7md.out: NSTEP = 61340 TIME(PS) = 1061.340 TEMP(K) = 391.71
> PRESS = 0.0
> 7md.out: NSTEP = 61350 TIME(PS) = 1061.350 TEMP(K) = 387.40
> PRESS = 0.0
> 7md.out: NSTEP = 61360 TIME(PS) = 1061.360 TEMP(K) = 378.63
> PRESS = 0.0
> 7md.out: NSTEP = 61370 TIME(PS) = 1061.370 TEMP(K) = 365.20
> PRESS = 0.0
> 7md.out: NSTEP = 61380 TIME(PS) = 1061.380 TEMP(K) = 350.35
> PRESS = 0.0
> 7md.out: NSTEP = 61390 TIME(PS) = 1061.390 TEMP(K) = 337.03
> PRESS = 0.0
> 7md.out: NSTEP = 61400 TIME(PS) = 1061.400 TEMP(K) = 325.53
> PRESS = 0.0
> 7md.out: NSTEP = 61410 TIME(PS) = 1061.410 TEMP(K) = 319.78
> PRESS = 0.0
> 7md.out: NSTEP = 61420 TIME(PS) = 1061.420 TEMP(K) = 317.17
> PRESS = 0.0
> 7md.out: NSTEP = 61430 TIME(PS) = 1061.430 TEMP(K) = 323.28
> PRESS = 0.0
> 7md.out: NSTEP = 61440 TIME(PS) = 1061.440 TEMP(K) = 332.56
> PRESS = 0.0
> 7md.out: NSTEP = 61450 TIME(PS) = 1061.450 TEMP(K) = 345.61
> PRESS = 0.0
> 7md.out: NSTEP = 61460 TIME(PS) = 1061.460 TEMP(K) = NaN
> PRESS = 0.0
> 7md.out: NSTEP = 61470 TIME(PS) = 1061.470 TEMP(K) = NaN
> PRESS = 0.0
> 7md.out: NSTEP = 61480 TIME(PS) = 1061.480 TEMP(K) = NaN
> PRESS = 0.0
>
> The corresponding 7md.out log file has this when it crashes:
>
> 7md.out: NSTEP = 61460 TIME(PS) = 1061.460 TEMP(K) = NaN
> PRESS = 0.0
> 7md.out: Etot = NaN EKtot = NaN EPtot =
> -73908.6279
> 7md.out: NSTEP = 61470 TIME(PS) = 1061.470 TEMP(K) = NaN
> PRESS = 0.0
> 7md.out: Etot = NaN EKtot = NaN EPtot =
> **************
> 7md.out: NSTEP = 61480 TIME(PS) = 1061.480 TEMP(K) = NaN
> PRESS = 0.0
> 7md.out: Etot = NaN EKtot = NaN EPtot =
> **************
>
>
> Can you please guide me about the source of this error that TEMP and EKtot
> becomes NaN? I am following the general protocol described here
> https://ambermd.org/tutorials/basic/tutorial13/index.php
>
> Thank you for your time and help.
>
> Sincerely,
> PN
>
>
> On Wed, Feb 8, 2023 at 6:49 PM Prithviraj Nandigrami <
> prithviraj.nandigrami.gmail.com> wrote:
>
>> Thank you, Prof. Simmerling for your helpful response. Following your
>> suggestions, I was able to move past the minimization and heating stages
>> and now continue the subsequent steps of the relaxation protocol you
>> pointed out.
>>
>> One other thing I wanted to confirm. For nucleic acid, when I choose the
>> backbone to apply a restraint during relaxation, the following selection
>> mask should work?
>>
>> restraintmask=".OP1,OP2,OP3,P,C1',C2',C3',C4',C5',O2',O3',O4',O5'"
>>
>> The minimization in 5min.in looks like the following:
>>
>> Minimization of everything excluding backbone
>> &cntrl
>> imin = 1, maxcyc = 1000, ntmin = 2,
>> ncyc = 30, ntx = 1,
>> ntwe = 0, ntwr = 500, ntpr = 50,
>> ntc = 2, ntf = 2, ntb = 1, ntp = 0,
>> cut = 8.0,
>> ntr=1,
>> restraintmask=".OP1,OP2,OP3,P,C1',C2',C3',C4',C5',O2',O3',O4',O5'",
>> restraint_wt=5.0,
>> ioutfm=1, ntxo=2,
>> /
>>
>>
>> I was able to run until 4md.in following the protocol you kindly
>> referred me to. However, in the subsequent minimization, it fails with the
>> following error in the 5min.out file:
>>
>>
>> Minimization of everything excluding backbone
>>
>> &cntrl
>>
>> imin = 1, maxcyc = 1000, ntmin = 2,
>>
>> ncyc = 30, ntx = 1,
>>
>> ntwe = 0, ntwr = 500, ntpr = 50,
>>
>> ntc = 2, ntf = 2, ntb = 1, ntp = 0,
>>
>> cut = 8.0,
>>
>> ntr=1,
>> restraintmask=".OP1,OP2,OP3,P,C1',C2',C3',C4',C5',O2',O3',O4',O5'", re
>> ioutfm=1, ntxo=2,
>>
>> /
>>
>>
>>
>>
>>
>> Note: ig = -1. Setting random seed to 373100 based on wallclock time in
>> microseconds.
>> | irandom = 1, using AMBER's internal random number generator (default).
>>
>>
>> | Conditional Compilation Defines Used:
>> | PUBFFT
>> | BINTRAJ
>> | EMIL
>>
>> | Largest sphere to fit in unit cell has radius = 42.443
>>
>> | New format PARM file being parsed.
>> | Version = 1.000 Date = 02/06/23 Time = 17:52:04
>>
>> | Note: 1-4 EEL scale factors are being read from the topology file.
>>
>> | Note: 1-4 VDW scale factors are being read from the topology file.
>> | Duplicated 0 dihedrals
>>
>> | Duplicated 0 dihedrals
>>
>>
>> --------------------------------------------------------------------------------
>> 1. RESOURCE USE:
>>
>> --------------------------------------------------------------------------------
>>
>> getting box info from netcdf restart file
>> NATOM = 86021 NTYPES = 18 NBONH = 84544 MBONA = 1576
>> NTHETH = 1599 MTHETA = 2461 NPHIH = 3564 MPHIA = 4813
>> NHPARM = 0 NPARM = 0 NNB = 123624 NRES = 28078
>> NBONA = 1576 NTHETA = 2461 NPHIA = 4813 NUMBND = 50
>> NUMANG = 98 NPTRA = 62 NATYP = 31 NPHB = 1
>> IFBOX = 2 NMXRS = 36 IFCAP = 0 NEXTRA = 0
>> NCOPY = 0
>>
>> | Coordinate Index Table dimensions: 20 20 20
>> | Direct force subcell size = 5.1982 5.1982 5.1982
>>
>> BOX TYPE: TRUNCATED OCTAHEDRON
>>
>>
>> --------------------------------------------------------------------------------
>> 2. CONTROL DATA FOR THE RUN
>>
>> --------------------------------------------------------------------------------
>>
>> default_name
>>
>>
>> General flags:
>> imin = 1, nmropt = 0
>>
>> Nature and format of input:
>> ntx = 1, irest = 0, ntrx = 1
>>
>> Nature and format of output:
>> ntxo = 2, ntpr = 50, ntrx = 1, ntwr =
>> 500
>> iwrap = 0, ntwx = 0, ntwv = 0, ntwe =
>> 0
>> ioutfm = 1, ntwprt = 0, idecomp = 0, rbornstat=
>> 0
>>
>> Potential function:
>> ntf = 2, ntb = 1, igb = 0, nsnb =
>> 25
>> ipol = 0, gbsa = 0, iesp = 0
>> dielc = 1.00000, cut = 8.00000, intdiel = 1.00000
>>
>> Frozen or restrained atoms:
>> ibelly = 0, ntr = 1
>> restraint_wt = 5.00000
>>
>> Energy minimization:
>> maxcyc = 1000, ncyc = 30, ntmin = 1
>> dx0 = 0.01000, drms = 0.00010
>>
>> SHAKE:
>> ntc = 2, jfastw = 0
>> tol = 0.00001
>>
>> | Intermolecular bonds treatment:
>> | no_intermolecular_bonds = 1
>>
>> | Energy averages sample interval:
>> | ene_avg_sampling = 1
>>
>> Ewald parameters:
>> verbose = 0, ew_type = 0, nbflag = 1, use_pme =
>> 1
>> vdwmeth = 1, eedmeth = 1, netfrc = 0
>> Box X = 103.964 Box Y = 103.964 Box Z = 103.964
>> Alpha = 109.471 Beta = 109.471 Gamma = 109.471
>> NFFT1 = 108 NFFT2 = 108 NFFT3 = 108
>> Cutoff= 8.000 Tol =0.100E-04
>> Ewald Coefficient = 0.34864
>> Interpolation order = 4
>>
>> LOADING THE CONSTRAINED ATOMS AS GROUPS
>>
>>
>> 5. REFERENCE ATOM COORDINATES
>>
>> default_name
>>
>> Mask .OP1,OP2,OP3,P,C1',C2',C3',C4',C5',O2',O3',O4',O5'; matches
>> 781 atoms
>>
>>
>> --------------------------------------------------------------------------------
>> 3. ATOMIC COORDINATES AND VELOCITIES
>>
>> --------------------------------------------------------------------------------
>>
>> default_name
>>
>> begin time read from input coords = 4000.000 ps
>>
>>
>> Number of triangulated 3-point waters found: 27948
>>
>> Sum of charges from parm topology file = -0.00000005
>> Forcing neutrality...
>>
>> | Dynamic Memory, Types Used:
>> | Reals 2361857
>> | Integers 2446997
>>
>> | Nonbonded Pairs Initial Allocation: 14367657
>>
>>
>> --------------------------------------------------------------------------------
>> 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.2738E-14 at 2.422500
>> | CHECK d/dx switch(x): max rel err = 0.8332E-11 at 2.782960
>> ---------------------------------------------------
>> |---------------------------------------------------
>> | APPROXIMATING direct energy using CUBIC SPLINE INTERPOLATION
>> | with 50.0 points per unit in tabled values
>> | Relative Error Limit not exceeded for r .gt. 2.47
>> | APPROXIMATING direct force using CUBIC SPLINE INTERPOLATION
>> | with 50.0 points per unit in tabled values
>> | Relative Error Limit not exceeded for r .gt. 2.89
>> |---------------------------------------------------
>>
>>
>> NSTEP ENERGY RMS GMAX NAME NUMBER
>> 1 -2.9042E+05 1.7342E+01 1.0165E+02 C6 2061
>>
>> BOND = 583.7031 ANGLE = 1084.0410 DIHED =
>> 1636.3989
>> VDWAALS = 38858.5291 EEL = -324802.9124 HBOND =
>> 0.0000
>> 1-4 VDW = 582.9706 1-4 EEL = -8366.0474 RESTRAINT =
>> 0.0000
>>
>> .... RESTARTED DUE TO LINMIN FAILURE ...
>>
>>
>> NSTEP ENERGY RMS GMAX NAME NUMBER
>> 50 -3.0204E+05 1.6344E+01 5.4097E+01 O 52158
>>
>> BOND = 139.4672 ANGLE = 514.9922 DIHED =
>> 1512.5095
>> VDWAALS = 30012.0788 EEL = -326465.3379 HBOND =
>> 0.0000
>> 1-4 VDW = 546.7331 1-4 EEL = -8313.5239 RESTRAINT =
>> 15.5416
>> EAMBER = -302053.0811
>>
>> .... RESTARTED DUE TO LINMIN FAILURE ...
>>
>> .... RESTARTED DUE TO LINMIN FAILURE ...
>>
>>
>> NSTEP ENERGY RMS GMAX NAME NUMBER
>> 100 -3.0211E+05 1.6309E+01 5.3292E+01 O 52158
>>
>> BOND = 130.3847 ANGLE = 480.4963 DIHED =
>> 1503.1108
>> VDWAALS = 29584.1699 EEL = -326052.3685 HBOND =
>> 0.0000
>> 1-4 VDW = 544.0760 1-4 EEL = -8315.4291 RESTRAINT =
>> 14.5677
>> EAMBER = -302125.5599
>>
>> .... RESTARTED DUE TO LINMIN FAILURE ...
>>
>> .... RESTARTED DUE TO LINMIN FAILURE ...
>>
>>
>> FINAL RESULTS
>>
>>
>>
>> NSTEP ENERGY RMS GMAX NAME NUMBER
>> 139 -3.0211E+05 1.6309E+01 5.3292E+01 O 52158
>>
>> BOND = 130.3846 ANGLE = 480.4963 DIHED =
>> 1503.1108
>> VDWAALS = 29584.1683 EEL = -326052.3666 HBOND =
>> 0.0000
>> 1-4 VDW = 544.0760 1-4 EEL = -8315.4291 RESTRAINT =
>> 14.5677
>> EAMBER = -302125.5597
>>
>> ***** REPEATED LINMIN FAILURE *****
>>
>>
>> Is there anything different I need to do at this stage because it
>> encounters a LINMIN error?
>>
>> Thank you for your time and help.
>>
>> Sincerely,
>> PN
>>
>>
>> On Wed, Feb 8, 2023 at 3:53 PM Carlos Simmerling <
>> carlos.simmerling.gmail.com> wrote:
>>
>>> I'm not sure, nothing seems obviously wrong. there are a few things you
>>> might try:
>>>
>>> - I don't usually put restraints on H atoms
>>> - you could try fewer minimization steps, I usually only do 200 or
>>> so, with ntmin=2.
>>> - try printing more often (ntpr = 10, or even 1) when there are
>>> errors, it can help you see what's going on.
>>> - look at your structure near the atoms with high gradient (gmax, P
>>> 559 in your output).
>>>
>>>
>>> On Wed, Feb 8, 2023 at 2:54 PM Prithviraj Nandigrami <
>>> prithviraj.nandigrami.gmail.com> wrote:
>>>
>>>> Thank you, Prof. Simmerling for pointing out the more modern relaxation
>>>> protocol.
>>>>
>>>> I used the following input file for the first minimization for solvent:
>>>>
>>>> minimization of solvent
>>>> &cntrl
>>>> imin = 1, maxcyc = 1000,
>>>> ncyc = 20, ntx = 1,
>>>> ntwe = 0, ntwr = 500, ntpr = 50,
>>>> ntc = 1, ntf = 2, ntb = 1, ntp = 0,
>>>> cut = 10.0,
>>>> ntr=1, restraintmask = ':1-65',
>>>> restraint_wt = 10.,
>>>> ioutfm=1, ntxo=2,
>>>> &end
>>>>
>>>> However, it soon crashes with the following error message:
>>>>
>>>> Program received signal SIGSEGV: Segmentation fault - invalid memory
>>>> reference.
>>>>
>>>> Backtrace for this error:
>>>> #0 0x7f0d1dc83e3a
>>>> #1 0x7f0d1dc83013
>>>> #2 0x7f0d1d47d3ff
>>>> #3 0x4bbdd2
>>>> #4 0x4b1e13
>>>> #5 0x4e196a
>>>> #6 0x50fff8
>>>> #7 0x42199c
>>>> #8 0x7f0d1d469554
>>>> #9 0x4219cc
>>>> #10 0xffffffffffffffff
>>>> Segmentation fault (core dumped)
>>>>
>>>> When I look at the 1min.out, I see the following:
>>>>
>>>> Here is the input file:
>>>>
>>>> minimization of solvent
>>>>
>>>> &cntrl
>>>>
>>>> imin = 1, maxcyc = 1000,
>>>>
>>>> ncyc = 20, ntx = 1,
>>>>
>>>> ntwe = 0, ntwr = 500, ntpr = 50,
>>>>
>>>> ntc = 1, ntf = 2, ntb = 1, ntp = 0,
>>>>
>>>> cut = 10.0,
>>>>
>>>> ntr=1, restraintmask = ':1-65',
>>>>
>>>> restraint_wt = 10.,
>>>>
>>>> ioutfm=1, ntxo=2,
>>>>
>>>> &end
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> Note: ig = -1. Setting random seed to 482884 based on wallclock time
>>>> in
>>>> microseconds.
>>>> | irandom = 1, using AMBER's internal random number generator (default).
>>>>
>>>>
>>>> | Conditional Compilation Defines Used:
>>>> | PUBFFT
>>>> | BINTRAJ
>>>> | EMIL
>>>>
>>>> | Largest sphere to fit in unit cell has radius = 43.564
>>>>
>>>> | New format PARM file being parsed.
>>>> | Version = 1.000 Date = 02/06/23 Time = 17:52:04
>>>>
>>>> | Note: 1-4 EEL scale factors are being read from the topology file.
>>>>
>>>> | Note: 1-4 VDW scale factors are being read from the topology file.
>>>> | Duplicated 0 dihedrals
>>>>
>>>> | Duplicated 0 dihedrals
>>>>
>>>>
>>>> --------------------------------------------------------------------------------
>>>> 1. RESOURCE USE:
>>>>
>>>> --------------------------------------------------------------------------------
>>>>
>>>> getting new box info from bottom of inpcrd
>>>> NATOM = 86021 NTYPES = 18 NBONH = 84544 MBONA = 1576
>>>> NTHETH = 1599 MTHETA = 2461 NPHIH = 3564 MPHIA = 4813
>>>> NHPARM = 0 NPARM = 0 NNB = 123624 NRES = 28078
>>>> NBONA = 1576 NTHETA = 2461 NPHIA = 4813 NUMBND = 50
>>>> NUMANG = 98 NPTRA = 62 NATYP = 31 NPHB = 1
>>>> IFBOX = 2 NMXRS = 36 IFCAP = 0 NEXTRA = 0
>>>> NCOPY = 0
>>>>
>>>> | Coordinate Index Table dimensions: 17 17 17
>>>> | Direct force subcell size = 6.2771 6.2771 6.2771
>>>>
>>>> BOX TYPE: TRUNCATED OCTAHEDRON
>>>>
>>>>
>>>> --------------------------------------------------------------------------------
>>>> 2. CONTROL DATA FOR THE RUN
>>>>
>>>> --------------------------------------------------------------------------------
>>>>
>>>> default_name
>>>>
>>>>
>>>> General flags:
>>>> imin = 1, nmropt = 0
>>>>
>>>> Nature and format of input:
>>>> ntx = 1, irest = 0, ntrx = 1
>>>>
>>>> Nature and format of output:
>>>> ntxo = 2, ntpr = 50, ntrx = 1, ntwr =
>>>> 500
>>>> iwrap = 0, ntwx = 0, ntwv = 0, ntwe =
>>>> 0
>>>> ioutfm = 1, ntwprt = 0, idecomp = 0,
>>>> rbornstat= 0
>>>>
>>>> Potential function:
>>>> ntf = 2, ntb = 1, igb = 0, nsnb =
>>>> 25
>>>> ipol = 0, gbsa = 0, iesp = 0
>>>> dielc = 1.00000, cut = 10.00000, intdiel = 1.00000
>>>>
>>>> Frozen or restrained atoms:
>>>> ibelly = 0, ntr = 1
>>>> restraint_wt = 10.00000
>>>>
>>>> Energy minimization:
>>>> maxcyc = 1000, ncyc = 20, ntmin = 1
>>>> dx0 = 0.01000, drms = 0.00010
>>>>
>>>> | Intermolecular bonds treatment:
>>>> | no_intermolecular_bonds = 1
>>>>
>>>> | Energy averages sample interval:
>>>> | ene_avg_sampling = 1
>>>>
>>>> Ewald parameters:
>>>> verbose = 0, ew_type = 0, nbflag = 1, use_pme =
>>>> 1
>>>> vdwmeth = 1, eedmeth = 1, netfrc = 0
>>>> Box X = 106.711 Box Y = 106.711 Box Z = 106.711
>>>> Alpha = 109.471 Beta = 109.471 Gamma = 109.471
>>>> NFFT1 = 108 NFFT2 = 108 NFFT3 = 108
>>>> Cutoff= 10.000 Tol =0.100E-04
>>>> Ewald Coefficient = 0.27511
>>>> Interpolation order = 4
>>>>
>>>> LOADING THE CONSTRAINED ATOMS AS GROUPS
>>>>
>>>>
>>>> 5. REFERENCE ATOM COORDINATES
>>>>
>>>> default_name
>>>>
>>>> Mask :1-65; matches 2112 atoms
>>>>
>>>>
>>>> --------------------------------------------------------------------------------
>>>> 3. ATOMIC COORDINATES AND VELOCITIES
>>>>
>>>> --------------------------------------------------------------------------------
>>>>
>>>> default_name
>>>>
>>>> begin time read from input coords = 0.000 ps
>>>>
>>>>
>>>> Number of triangulated 3-point waters found: 27948
>>>>
>>>> Sum of charges from parm topology file = -0.00000005
>>>> Forcing neutrality...
>>>>
>>>> | Dynamic Memory, Types Used:
>>>> | Reals 2381613
>>>> | Integers 2445597
>>>>
>>>> | Nonbonded Pairs Initial Allocation: 26019201
>>>>
>>>>
>>>> --------------------------------------------------------------------------------
>>>> 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.2738E-14 at 2.422500
>>>> | CHECK d/dx switch(x): max rel err = 0.8314E-11 at 2.736960
>>>> ---------------------------------------------------
>>>> |---------------------------------------------------
>>>> | APPROXIMATING direct energy using CUBIC SPLINE INTERPOLATION
>>>> | with 50.0 points per unit in tabled values
>>>> | Relative Error Limit not exceeded for r .gt. 2.33
>>>> | APPROXIMATING direct force using CUBIC SPLINE INTERPOLATION
>>>> | with 50.0 points per unit in tabled values
>>>> | Relative Error Limit not exceeded for r .gt. 2.80
>>>> |---------------------------------------------------
>>>>
>>>>
>>>> NSTEP ENERGY RMS GMAX NAME
>>>> NUMBER
>>>> 1 -2.4046E+05 1.7331E+02 3.5757E+04 P 559
>>>>
>>>> BOND = 155.7912 ANGLE = 941.4734 DIHED =
>>>> 1512.3716
>>>> VDWAALS = 46273.0244 EEL = -282423.9933 HBOND =
>>>> 0.0000
>>>> 1-4 VDW = 823.9232 1-4 EEL = -7741.5257 RESTRAINT =
>>>> 0.0000
>>>>
>>>> ----------------------------------------------------------
>>>>
>>>> What am I doing incorrectly?
>>>>
>>>> On Wed, Feb 8, 2023 at 11:10 AM Carlos Simmerling via AMBER <
>>>> amber.ambermd.org> wrote:
>>>>
>>>>> I'm not sure about later steps, but your heating input has errors. You
>>>>> set
>>>>> nstlim=1000, but then for the thermostat changes you have this:
>>>>>
>>>>> &wt type='TEMP0',istep1=0,
>>>>>
>>>>> istep2=500, value1=100.0,value2=300.0,
>>>>>
>>>>> &end
>>>>>
>>>>> &wt type='TEMP0',istep1=5001,
>>>>>
>>>>> istep2=10000, value1=300.0,value2=300.0,
>>>>>
>>>>> &end
>>>>>
>>>>> &wt type='END', &end
>>>>>
>>>>>
>>>>> notice that you heat very quickly (500 MD steps) to 300K, and then
>>>>> define
>>>>> it again for step 5001 and on, which does not exist in this simulation.
>>>>> if you look at the mdout file, you can see that at step 500 it has not
>>>>> yet
>>>>> reached 300K. This may be from using ntt=1, which is not the thermostat
>>>>> that people normally use (I think ntt=3 is much more common). You
>>>>> could
>>>>> also modify the coupling constant (tautp), but I think you're better
>>>>> off
>>>>> using ntt=3. The temperature continues to slowly rise, but has not
>>>>> reached
>>>>> your target at the end, so you should not continue to the next step
>>>>> yet.
>>>>>
>>>>> overall you are using a relatively simple relaxation protocol. I
>>>>> suggest
>>>>> looking at the Amber web page tutorials to see something that is more
>>>>> modern. it isn't for RNA, but you should be able to adapt the general
>>>>> ideas.
>>>>> https://ambermd.org/tutorials/basic/tutorial13/index.php
>>>>>
>>>>> On Wed, Feb 8, 2023 at 10:47 AM Prithviraj Nandigrami <
>>>>> prithviraj.nandigrami.gmail.com> wrote:
>>>>>
>>>>> > Thank you, Prof. Simmerling for your helpful response. I apologize
>>>>> for
>>>>> > sending the wrong output files. Indeed, it does look like I was able
>>>>> to
>>>>> > solve the issue of low density by modifying how I built the waterbox
>>>>> > following your suggestion in an earlier email.
>>>>> >
>>>>> > However, I notice now that there is a new issue. In the equilibration
>>>>> > stage, there is an intermittent error of exceeding 'vlimit' as below:
>>>>> >
>>>>> > eq.mdout: t = 0.00000, dt = 0.00200, vlimit =
>>>>> 20.00000
>>>>> > eq.mdout:vlimit exceeded for step 67130; vmax = 41.0592
>>>>> > eq.mdout:vlimit exceeded for step 67131; vmax = 31.1308
>>>>> > eq.mdout:vlimit exceeded for step 67132; vmax = 29.9301
>>>>> > eq.mdout:vlimit exceeded for step 67133; vmax = 20.7105
>>>>> > eq.mdout:vlimit exceeded for step 67135; vmax = 26.6603
>>>>> > eq.mdout:vlimit exceeded for step 67136; vmax = 38.0706
>>>>> > eq.mdout:vlimit exceeded for step 67145; vmax = 21.1061
>>>>> > eq.mdout:vlimit exceeded for step 67206; vmax = 21.1132
>>>>> > eq.mdout:vlimit exceeded for step 67217; vmax = 20.3974
>>>>> > eq.mdout:vlimit exceeded for step 67276; vmax = 22.7971
>>>>> > eq.mdout:vlimit exceeded for step 67298; vmax = 21.0155
>>>>> > eq.mdout:vlimit exceeded for step 67302; vmax = 20.4616
>>>>> > eq.mdout:vlimit exceeded for step 67324; vmax = 20.8580
>>>>> > eq.mdout:vlimit exceeded for step 110108; vmax = 35.5854
>>>>> > eq.mdout:vlimit exceeded for step 110109; vmax = 22.3098
>>>>> > eq.mdout:vlimit exceeded for step 110110; vmax = 23.9470
>>>>> > eq.mdout:vlimit exceeded for step 110113; vmax = 21.0085
>>>>> > eq.mdout:vlimit exceeded for step 110114; vmax = 21.2495
>>>>> > eq.mdout:vlimit exceeded for step 110116; vmax = 21.5208
>>>>> > eq.mdout:vlimit exceeded for step 110124; vmax = 221.6161
>>>>> > eq.mdout:vlimit exceeded for step 110125; vmax = 25.9272
>>>>> > eq.mdout:vlimit exceeded for step 110126; vmax = 20.7325
>>>>> > eq.mdout:vlimit exceeded for step 110127; vmax = 24.5213
>>>>> > eq.mdout:vlimit exceeded for step 110129; vmax = 21.0525
>>>>> > eq.mdout:vlimit exceeded for step 110131; vmax = 23.6221
>>>>> > eq.mdout:vlimit exceeded for step 110133; vmax = 28.7984
>>>>> >
>>>>> > I do not see this error in the heating stage. However, the subsequent
>>>>> > production run crashes. I attached the input files I used and the
>>>>> > corresponding output files. In the md_1.mdout file, after a certain
>>>>> number
>>>>> > of steps (at TIME = 237.54 ps), I notice that TEMP(K) = NaN, and
>>>>> PRESS =
>>>>> > ******
>>>>> > What is odd is that the value of PRESS as it is printed out
>>>>> fluctuates
>>>>> > between positive and negative values. Does this indicate a problem
>>>>> with any
>>>>> > input file parameter?
>>>>> >
>>>>> > Thank you for your time and continued help.
>>>>> >
>>>>> > Sincerely,
>>>>> > PN
>>>>> >
>>>>> > eq.mdout
>>>>> > <
>>>>> https://drive.google.com/file/d/10GDaOG7N-iSBTinwoOvQLdHF-4k7K2ce/view?usp=drive_web
>>>>> >
>>>>> > heat.mdout
>>>>> > <
>>>>> https://drive.google.com/file/d/1A9HrHr2vxWz9lI5Y6QHQbTM6GOdQia-S/view?usp=drive_web
>>>>> >
>>>>> > md_1.mdout
>>>>> > <
>>>>> https://drive.google.com/file/d/1RrbCeCc9nq9MvpPRaD6okvP_olOGAloK/view?usp=drive_web
>>>>> >
>>>>> > mini.mdout
>>>>> > <
>>>>> https://drive.google.com/file/d/1xEdo65b_kSH6J3Ielmr5hNIlR4g1NF9H/view?usp=drive_web
>>>>> >
>>>>> >
>>>>> >
>>>>> >
>>>>> > On Wed, Feb 8, 2023 at 7:21 AM Carlos Simmerling via AMBER <
>>>>> > amber.ambermd.org> wrote:
>>>>> >
>>>>> >> the "step4.1_equilibration" output here doesn't match what you sent
>>>>> >> before.
>>>>> >> this one goes to a density of 1 rather quickly, which Prof Case
>>>>> mentioned
>>>>> >> was a problem in the previous one. Maybe you need to send updated
>>>>> >> information about the problem observed with these specific outputs
>>>>> - I
>>>>> >> didn't see any error in the files that you attached.
>>>>> >>
>>>>> >> On Tue, Feb 7, 2023 at 10:37 AM Prithviraj Nandigrami <
>>>>> >> prithviraj.nandigrami.gmail.com> wrote:
>>>>> >>
>>>>> >> > Thank you, Dr. Case and Dr. Simmerling for your helpful
>>>>> suggestions. The
>>>>> >> > input pdb file was processed through 'pdb4amber' prior to system
>>>>> >> > preparation using tleap.
>>>>> >> >
>>>>> >> > Here are the commands I used to prepare and solvate the system in
>>>>> tleap:
>>>>> >> >
>>>>> >> > parm10 = loadamberparams parm10.dat
>>>>> >> > loadOff terminal_monophosphate.lib
>>>>> >> > source leaprc.water.tip3p
>>>>> >> > loadAmberParams frcmod.tip3p
>>>>> >> > loadamberparams M6A.frcmod
>>>>> >> > M6A = loadMol2 M6A.mol2
>>>>> >> > addions mol Na+ 0
>>>>> >> > solvateOct mol TIP3PBOX 12.0 0.75
>>>>> >> > saveAmberParm mol 2L1F.GAACA.M6A.parm7 2L1F.GAACA.M6A.rst7
>>>>> >> >
>>>>> >> > I also attached the pdb4amber output file (before tleap) and the
>>>>> output
>>>>> >> > prepared system after tleap, as well as the output log file for
>>>>> tleap.
>>>>> >> > Also, attached is the AMBER minimization and equilibration log
>>>>> files.
>>>>> >> >
>>>>> >> > I would appreciate any help regarding this.
>>>>> >> >
>>>>> >> > Sincerely,
>>>>> >> > PN
>>>>> >> >
>>>>> >> >
>>>>> >> >
>>>>> >> >
>>>>> >> >
>>>>> >> >
>>>>> >> > On Tue, Feb 7, 2023 at 8:26 AM David A Case <
>>>>> david.case.rutgers.edu>
>>>>> >> > wrote:
>>>>> >> >
>>>>> >> >> On Mon, Feb 06, 2023, Prithviraj Nandigrami via AMBER wrote:
>>>>> >> >> >This is what the density looks like in the mdout files:
>>>>> >> >> >
>>>>> >> >> >*Equilibration *
>>>>> >> >> >
>>>>> >> >> > 0.8121
>>>>> >> >> > 0.7651
>>>>> >> >> > 0.7322
>>>>> >> >> > 0.7139
>>>>> >> >> > 0.7063
>>>>> >> >> ....
>>>>> >> >> >*Production*
>>>>> >> >> >
>>>>> >> >> > 0.7694
>>>>> >> >> > 0.7699
>>>>> >> >> > 0.7722
>>>>> >> >> > 0.7744
>>>>> >> >> > 0.7767
>>>>> >> >>
>>>>> >> >> ...
>>>>> >> >>
>>>>> >> >> These are *very* low densities, certainly for a simulation of a
>>>>> >> >> macromolecule in water. The fact that they don't quickly
>>>>> approach a
>>>>> >> value
>>>>> >> >> near 1.0 (since most of a typical system is water) suggests some
>>>>> >> problem
>>>>> >> >> with the setup of the system. I don't have enough info now to
>>>>> grok
>>>>> >> >> what the most likely problem is there.
>>>>> >> >>
>>>>> >> >> ....dac
>>>>> >> >>
>>>>> >> >>
>>>>> >> _______________________________________________
>>>>> >> AMBER mailing list
>>>>> >> AMBER.ambermd.org
>>>>> >> http://lists.ambermd.org/mailman/listinfo/amber
>>>>> >>
>>>>> >
>>>>> _______________________________________________
>>>>> AMBER mailing list
>>>>> AMBER.ambermd.org
>>>>> http://lists.ambermd.org/mailman/listinfo/amber
>>>>>
>>>>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Feb 09 2023 - 08:00:02 PST
