that it something you will need to explore. my suspicion would be that the
reduced restraints allow more motion, which reveals a problem in the force
field. looking at the structure snapshots when the temperature is unstable
may help, as would examining the new force field parameters. When a
simulation is unstable with modified parameters, that is the most likely
source of the problem. it could also be that the coordinates are unreliable.
On Thu, Feb 9, 2023 at 10:46 AM Prithviraj Nandigrami <
prithviraj.nandigrami.gmail.com> wrote:
> Thank you, Prof. Simmerling for your reply.
>
> The input PDB file has a modified methylated nucleotide which was built
> using tleap steps prior. What I am confused about is that the previous step
> of the protocol runs fine, with the only difference between the previous
> step and the step where it crashes is that the value of restraint_wt is
> lower (10.0 vs 1.0).
>
> On Thu, Feb 9, 2023 at 10:36 AM Carlos Simmerling via AMBER <
> amber.ambermd.org> wrote:
>
>> 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
>>
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Feb 09 2023 - 08:30:05 PST
