Re: FATAL ERROR: Periodic cell has become too small for original patch grid!

From: Peter Freddolino (petefred_at_ks.uiuc.edu)
Date: Wed Jan 30 2008 - 15:57:53 CST

Hi Josh,
if I recall correctly your original protocol was:

-prepare system
-run a set of NVE simulations with constraints, with velocity rescaling
every some number of steps
-run an NPT simulation with gentler constraints

Can I verify that your protocol is now:
-prepare system
-run a set of NPT simulations with constraints, gradually heating the system
-run an NPT simulation with gentler constraints

Is this correct? How does the cell size look during that first simulation?
Peter

Joshua Adelman wrote:
> Hi Peter,
>
> I'm still running into the same problems. I also tried minimizing the
> system longer, since I was doing a rather small number of steps
> (originally 200 and now 1500), and have changed the piston parameters
> as you suggested. Also I added a line inside the loop that also set
> the langevin piston temperature to the same temp as the thermostat temp.
>
> I also just resolvated the system using autoionize (the first time I
> used meadionize), and I'm still getting the same error. I also tried
> just restraining the CA's and I'm still getting the same. Each thing
> I've tried, I've started from the minimization. When I get rid of the
> constraints, everything seems to work fine.
>
> I'm baffled . . .
>
> Josh
>
>
> On Jan 30, 2008, at 9:55 AM, Peter Freddolino wrote:
>
>> Hi Josh,
>> yes, although you can also wrap it in a tcl loop to make it easier:
>>
>> for {set i 15} {$i < 315} {incr i 15} {
>> set temperature $i
>> langevinTemp $temperature
>> run 1000
>> }
>>
>> (you could also add rescalevels statements if you're so inclined)
>>
>> Peter
>>
>> Joshua Adelman wrote:
>>> Hi Peter,
>>>
>>> Can you give me an example of how to implement the heating protocol
>>> using langevinTemp? Your thermostat would have something like:
>>> # Constant Temperature Control
>>> langevin on ;# do langevin dynamics
>>> langevinDamping 5 ;# damping coefficient (gamma) of 5/ps
>>> langevinTemp $temperature
>>> langevinHydrogen no ;# don't couple langevin bath to hydrogens
>>>
>>> Would you just tack the following on the end of your input script?
>>> set temperature 15
>>> langevinTemp $temperature
>>> run 1000
>>>
>>> set temperature 30
>>> langevinTemp $temperature
>>> run 1000
>>>
>>> set temperature 45
>>> langevinTemp $temperature
>>> run 1000
>>>
>>> etc.?
>>>
>>> Josh
>>>
>>> On Jan 30, 2008, at 9:34 AM, Peter Freddolino wrote:
>>>
>>>> Hi Josh,
>>>> actually, what I'd recommend instead is doing your initial heating
>>>> in NPT; *more* NV(E) equilibration is unlikely to help because it
>>>> is actually exacerbating the problem that because solvation
>>>> protocols are imperfect, we usually end up with a box that needs to
>>>> shrink somewhat; if you then feed the barostat a box that has
>>>> shrunk, there's a lot of vacuum and bad things start to happen.
>>>>
>>>> If you have constraints on your solute doing the heating in NPT is
>>>> unlikely to hurt the system; you can change the langevin dynamics
>>>> temperature with langevinTemp $NEWTEMP and follow the same heating
>>>> protocol as before. You may also want to increase your piston
>>>> period and piston decay by a factor of 2; 100/50 and 200/100 are
>>>> both reasonable values, but the latter averages over a longer
>>>> amount of time and thus is more likely to tolerate rapid initial
>>>> fluctuations.
>>>>
>>>> BTW, there's no need for it now, but it's generally helpful to
>>>> continue these discussions on-list so that they're searchable in
>>>> the future.
>>>> Best,
>>>> Peter
>>>>
>>>> Joshua Adelman wrote:
>>>>> Hi Peter,
>>>>>
>>>>> Good call. I printed out the DCD and XST files for every step and
>>>>> the system is definitely blowing up:
>>>>> # NAMD extended system trajectory file
>>>>> #$LABELS step a_x a_y a_z b_x b_y b_z c_x c_y c_z o_x o_y o_z s_x
>>>>> s_y s_z s_u s_v s_w
>>>>> 35000 119.6 0 0 0 87.2 0 0 0 112 0 0 0 0 0 0 0 0 0
>>>>> 35001 119.6 0 0 0 87.2 0 0 0 112 0 0 0 -0.000284665 -0.000284665
>>>>> -0.000284665 0 0 0
>>>>> 35002 119.397 0 0 0 87.0524 0 0 0 111.81 0 0 0 0.000385176
>>>>> 0.000385176 0.000385176 0 0 0
>>>>> 35003 119.397 0 0 0 87.0524 0 0 0 111.81 0 0 0 0.00202167
>>>>> 0.00202167 0.00202167 0 0 0
>>>>> 35004 120.763 0 0 0 88.0481 0 0 0 113.089 0 0 0 -0.00277211
>>>>> -0.00277211 -0.00277211 0 0 0
>>>>> 35005 120.763 0 0 0 88.0481 0 0 0 113.089 0 0 0 -0.014042
>>>>> -0.014042 -0.014042 0 0 0
>>>>> 35006 111.611 0 0 0 81.3755 0 0 0 104.519 0 0 0 0.0167054
>>>>> 0.0167054 0.0167054 0 0 0
>>>>> 35007 111.611 0 0 0 81.3755 0 0 0 104.519 0 0 0 0.0898706
>>>>> 0.0898706 0.0898706 0 0 0
>>>>> 35008 185.206 0 0 0 135.033 0 0 0 173.437 0 0 0 -0.324938
>>>>> -0.324938 -0.324938 0 0 0
>>>>> 35009 185.206 0 0 0 135.033 0 0 0 173.437 0 0 0 -1.23019 -1.23019
>>>>> -1.23019 0 0 0
>>>>>
>>>>> Although it starts to expand a little then shrinks a lot and then
>>>>> explodes, I'm guessing because it got compressed and is reacting
>>>>> to that. I never get any errors about atoms moving too fast though.
>>>>>
>>>>> Would the solution be to equilibrate longer in NVT using the
>>>>> restraints before turning on the barostat? Additionally do you
>>>>> think it would help to reduce the restraints on non-backbone atoms
>>>>> while still in NVT, before switching to NPT (and then I would
>>>>> relax the backbone restraints in the NPT ensemble)?
>>>>>
>>>>> Thanks again for your help.
>>>>>
>>>>> Josh
>>>>>
>>>>> On Jan 30, 2008, at 8:52 AM, Peter Freddolino wrote:
>>>>>
>>>>>> Hi Josh,
>>>>>> you may want to check how immediate 'immediately' really is. To
>>>>>> do this, try the same NPT run as before, but set your restart,
>>>>>> output, and dcd frequencies. See if it's really crashing on step
>>>>>> 0, or if the periodic cell rapidly shrinks before the crash (have
>>>>>> a look at your xst file). I've seen similar behavior in the past
>>>>>> when the water box contracts too much, too quickly; you may also
>>>>>> want to look at your velocity rescaling trajectories and see
>>>>>> whether there are any voids opening up during that portion of the
>>>>>> simulation, which would leave the barostat with a lot to fix very
>>>>>> quickly once you turn on pressure control.
>>>>>> Best,
>>>>>> Peter
>>>>>>
>>>>>> Joshua Adelman wrote:
>>>>>>> I am try to equilibrate a system, by first using harmonic
>>>>>>> constraints (restraints) on the protein atoms. I run a short
>>>>>>> simulation in which I hold the protein atoms and let the water
>>>>>>> and ions go as I slowly heat the system using progressive
>>>>>>> temperature rescaling without any thermo- or barostats. Once I
>>>>>>> reach 300K, I try to run another simulation where I relax the
>>>>>>> restraints on the non-backbone atoms of the protein in the NPT
>>>>>>> ensemble, but the simulation immediately crashes out giving the
>>>>>>> error:
>>>>>>>
>>>>>>> FATAL ERROR: Periodic cell has become too small for original
>>>>>>> patch grid!
>>>>>>> Possible solutions are to restart from a recent checkpoint,
>>>>>>> increase margin, or disable useFlexibleCell for liquid simulation.
>>>>>>>
>>>>>>> I tried increasing the margin to 2.0, but it does the same
>>>>>>> thing. useFlexibleCell is off already. If I turn off the
>>>>>>> barostat (NVT), the simulation runs fine and if I remove the
>>>>>>> constraints completely, NPT runs fine. I'm not sure how the
>>>>>>> barostat and the constraints are interacting to cause this
>>>>>>> problem. I'm including the input script below.
>>>>>>>
>>>>>>> Any suggestions would be appreciated.
>>>>>>>
>>>>>>> Josh
>>>>>>>
>>>>>>> #############################################################
>>>>>>> ## ADJUSTABLE PARAMETERS ##
>>>>>>> #############################################################
>>>>>>>
>>>>>>>
>>>>>>> structure ../2sub2.psf
>>>>>>> coordinates ../2sub2.pdb
>>>>>>> outputName 2s2-equil4
>>>>>>>
>>>>>>> set temperature 300.0
>>>>>>>
>>>>>>> # Continuing a job from the restart files
>>>>>>> set inputname 2s2-equil3
>>>>>>> binCoordinates $inputname.restart.coor
>>>>>>> binVelocities $inputname.restart.vel ;# remove the
>>>>>>> "temperature" entry if you use this!
>>>>>>> extendedSystem $inputname.restart.xsc
>>>>>>>
>>>>>>> set xscfile $inputname.restart.xsc
>>>>>>>
>>>>>>> #########################################################
>>>>>>> proc get_first_ts { xscfile } {
>>>>>>> set fd [open $xscfile r]
>>>>>>> gets $fd
>>>>>>> gets $fd
>>>>>>> gets $fd line
>>>>>>> set ts [lindex $line 0]
>>>>>>> close $fd
>>>>>>> return $ts
>>>>>>> }
>>>>>>>
>>>>>>> set firstts [get_first_ts $xscfile]
>>>>>>> ########################################################
>>>>>>>
>>>>>>> firsttimestep $firstts
>>>>>>>
>>>>>>>
>>>>>>> #############################################################
>>>>>>> ## SIMULATION PARAMETERS ##
>>>>>>> #############################################################
>>>>>>>
>>>>>>> # Input
>>>>>>> paraTypeCharmm on
>>>>>>> parameters ../par_all27_prot_na.prm
>>>>>>>
>>>>>>>
>>>>>>> # NOTE: Do not set the initial velocity temperature if you
>>>>>>> # have also specified a .vel restart file!
>>>>>>>
>>>>>>>
>>>>>>> # Periodic Boundary conditions
>>>>>>> # NOTE: Do not set the periodic cell basis if you have also
>>>>>>> # specified an .xsc restart file!
>>>>>>> #cellBasisVector1 119.5 0. 0.
>>>>>>> #cellBasisVector2 0. 87. 0.
>>>>>>> #cellBasisVector3 0. 0 111.9
>>>>>>> cellOrigin 0. 0. 0.
>>>>>>>
>>>>>>> wrapWater on
>>>>>>> wrapAll off
>>>>>>>
>>>>>>>
>>>>>>> # Force-Field Parameters
>>>>>>> exclude scaled1-4
>>>>>>> 1-4scaling 1.0
>>>>>>> cutoff 12.
>>>>>>> switching on
>>>>>>> switchdist 10.
>>>>>>> pairlistdist 13.5
>>>>>>> COMmotion yes
>>>>>>> margin 2.0
>>>>>>>
>>>>>>> # Integrator Parameters
>>>>>>> timestep 2.0 ;# 2fs/step
>>>>>>> rigidBonds all
>>>>>>> nonbondedFreq 1
>>>>>>> fullElectFrequency 2
>>>>>>> stepspercycle 10
>>>>>>>
>>>>>>> #PME (for full-system periodic electrostatics)
>>>>>>> PME yes
>>>>>>> PMEGridSizeX 125
>>>>>>> PMEGridSizeY 125
>>>>>>> PMEGridSizeZ 125
>>>>>>>
>>>>>>> # Restraints
>>>>>>> constraints on
>>>>>>> consref 1VYM_2sub2_restrained.pdb
>>>>>>> conskfile 1VYM_2sub2_restrained.pdb
>>>>>>> conskcol B
>>>>>>>
>>>>>>> # Constant Temperature Control
>>>>>>> langevin on ;# do langevin dynamics
>>>>>>> langevinDamping 5 ;# damping coefficient (gamma) of 5/ps
>>>>>>> langevinTemp $temperature
>>>>>>> langevinHydrogen no ;# don't couple langevin bath to
>>>>>>> hydrogens
>>>>>>>
>>>>>>>
>>>>>>> # Constant Pressure Control (variable volume)
>>>>>>> useGroupPressure yes ;# needed for 2fs steps
>>>>>>> useFlexibleCell no ;# no for water box, yes for membrane
>>>>>>> useConstantArea no ;# no for water box, yes for membrane
>>>>>>>
>>>>>>> langevinPiston on
>>>>>>> langevinPistonTarget 1.01325 ;# in bar -> 1 atm
>>>>>>> langevinPistonPeriod 100.
>>>>>>> langevinPistonDecay 50.
>>>>>>> langevinPistonTemp $temperature
>>>>>>>
>>>>>>> restartfreq 500 ;# 1000steps = every 1ps
>>>>>>> dcdfreq 500
>>>>>>> xstFreq 500
>>>>>>> outputEnergies 100
>>>>>>> outputPressure 100
>>>>>>>
>>>>>>> #protocol
>>>>>>> run 5000
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> ------------------------------------------------------------------------------------------------------
>>>>>>>
>>>>>>>

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:49:13 CST