Re: periodic boundary condition in XY, hard wall in Z - errors

From: Maxim Belkin (mbelkin_at_ks.uiuc.edu)
Date: Tue Mar 15 2016 - 15:50:18 CDT

Olya, you can do semi-hard wall with tcl-bc. All you have to do is

# in VMD:
set sel [atomselect top "all"]
set ch [open "somefile.txt" w]
puts $ch [$sel get serial]
close $ch
$sel delete

# in NAMD config file:
# // not tested, but should work ;) //

tclBC on
tclBCScript {

  set ch [open "somefile.txt" r]
  if {[gets $ch line] != -1} {
    set allAtomSerials [split $line]
  }

  set upperZlimit 50.0 ; # z < 50
  set lowerZlimit -50.0 ; # z > -50
  set padding 2.0
  set springConstant -0.1

  wrapmode input

  proc calcforces {step unique} {
    if {$step % 100 == 0 } {cleardrops;}
    while {[nextatom]} {
      set atomid [getid]
      foreach {x y z} [getcoord] { break; }
      if {$z < $::upperZlimit - 10.0} { dropatom; continue;}
      if {$z > $::lowerZlimit + 10.0} { dropatom; continue;}
      set force {0. 0. 0.}
      set fz 0.0
      if {$z > $::upperZlimit - $::padding} { set fz [expr {$::springConstant * ($z - $::upperZlimit + $::padding) }]; }
      if {$z < $::lowerZlimit + $::padding} { set fz [expr {$::springConstant * ($z - $::lowerZlimit - $::padding) }]; }
      lset force 2 $fz
      addforce $force
   
    }
  }
}
> On Mar 15, 2016, at 14:48, Olya Kravchenko <ovkrav_at_gmail.com> wrote:
>
> Maxim,
>
> thank you for replying to my questions. To conclude this thread I
> wanted to add, that after running a huge number of jobs in both (what
> I think is) NVT and NPT I noticed that --
>
> 1) When I run NVT with CellBasis3 commented out my box elongates in Z
> direction, so commenting it out is not equivalent to setting up a hard
> wall. I think if I run it very long it will just keep stretching to
> infinity or the limit that's allowed in namd. The size in X and Y
> holds well.
>
> 2) I don't have a margin problem when I do NVT, it only happens when I
> run NPT and only when I comment out CellBasis3. I think that is
> because there is nothing in Z direction, so unless I use different
> method for setting up a wall there, the data will make no sense.
>
> I am trying to figure out MSM now, but that would be a subject of
> another thread.
>
> Olga
>
>
> On Mon, Mar 14, 2016 at 12:24 PM, Maxim Belkin <mbelkin_at_ks.uiuc.edu> wrote:
>> Hi Olya,
>>
>> Yes, I meant langevinPiston. lagevin is for temperature control. You have to use langevin because you are using some mysterious tclForces that might heat up your system.
>>
>> Maxim
>>
>>
>>> On Mar 11, 2016, at 21:42, Olya Kravchenko <ovkrav_at_gmail.com> wrote:
>>>
>>> Maxim,
>>>
>>> thanks!
>>>
>>> Do you mean "langevin on" or "LangevinPiston on"? I thought "langevin
>>> on" controls temperature. If not, I am confused. I saw a couple of
>>> sample configuration files marked as NVT simulation where langevin is
>>> set to on... Are you saying that for NVT I need to comment it out as
>>> well?
>>>
>>> If you meant LangevinPiston, I just commented it out and ran another
>>> job and my simulation doesn't do anything crazy now. I also used exact
>>> measurements for cellbasis vectors.
>>>
>>> Olga
>>>
>>> On Fri, Mar 11, 2016 at 4:31 PM, Maxim Belkin <mbelkin_at_ks.uiuc.edu> wrote:
>>>> Olya,
>>>>
>>>> NAMD config script below has "langevin on", so this is NPT. Again, NPT simulations are used to equilibrate pressure, usually to 1 atm, which is also what you have set in your script. Vacuum has 0 atm. I have hard time imagining what langevinPiston should do in such a case.
>>>>
>>>> Here is a useful message on the subject:
>>>> http://www.ks.uiuc.edu/Research/namd/mailing_list/namd-l.2011-2012/4142.html
>>>>
>>>> Maxim
>>>>
>>>>> On Mar 11, 2016, at 14:57, Olya Kravchenko <ovkrav_at_gmail.com> wrote:
>>>>>
>>>>> Maxim,
>>>>>
>>>>> thank you for replying. I am trying to figure out MSM right now and
>>>>> see if it makes difference, I'll try setting exact values for
>>>>> cellbasisvectors too.
>>>>>
>>>>> As for the hard wall: I expected to see expansion in the beginning
>>>>> that does not change further, as it happens when I use PBC in 3D. What
>>>>> I saw is oscillating continuous expansion, i.e. a couple of times
>>>>> during simulation the whole box shrinks and then expands more than it
>>>>> did previous time; the change in proportion (elongation in z) is
>>>>> continuous.
>>>>>
>>>>> I ran simulations in NVT today just to see if it helps me keep the
>>>>> size of the box fixed and I still see the change in proportions and
>>>>> increase in Z, actually the whole box becomes much bigger, and
>>>>> proportions change as well, although not as much as it happens when I
>>>>> run NPT.
>>>>>
>>>>> Here is my configuration file for NVT:
>>>>>
>>>>> # Minimization and Equilibration of argon with PBC in XY
>>>>>
>>>>> structure ar.psf
>>>>> coordinates ar.pdb
>>>>> set outputname ar_force_out
>>>>>
>>>>> firsttimestep 0
>>>>>
>>>>> set temperature 310
>>>>>
>>>>> # Input
>>>>> paraTypeCharmm on
>>>>> parameters My_parameters.rtf
>>>>> temperature $temperature
>>>>>
>>>>>
>>>>> # Force-Field Parameters
>>>>> exclude scaled1-4
>>>>> 1-4scaling 1.0
>>>>> cutoff 12.0
>>>>> switching on
>>>>> switchdist 10.0
>>>>> pairlistdist 14.0
>>>>>
>>>>>
>>>>> # Integrator Parameters
>>>>> timestep 2.0
>>>>> rigidBonds none
>>>>> nonbondedFreq 1
>>>>> fullElectFrequency 2
>>>>> stepspercycle 10
>>>>>
>>>>>
>>>>> # Constant Temperature Control
>>>>> langevin on
>>>>> langevinDamping 5
>>>>> langevinTemp $temperature
>>>>> langevinHydrogen off
>>>>>
>>>>>
>>>>> # Periodic Boundary Conditions
>>>>> cellBasisVector1 62.0 0. 0.0
>>>>> cellBasisVector2 0.0 62.0 0.0
>>>>> #cellBasisVector3 0.0 0 0.0
>>>>> cellOrigin 30.0 30.0 32.5
>>>>>
>>>>> wrapAll on
>>>>>
>>>>> margin 30
>>>>>
>>>>> # Constant Pressure Control (variable volume)
>>>>> useGroupPressure no ;# needed for rigidBonds
>>>>> useFlexibleCell no
>>>>> useConstantArea no
>>>>>
>>>>> langevinPiston on
>>>>> langevinPistonTarget 1.01325 ;# in bar -> 1 atm
>>>>> langevinPistonPeriod 100.0
>>>>> langevinPistonDecay 50.0
>>>>> langevinPistonTemp $temperature
>>>>>
>>>>> tclForces on
>>>>> tclForcesScript my_force.tcl
>>>>>
>>>>> outputName $outputname
>>>>>
>>>>> restartfreq 500
>>>>> dcdfreq 250
>>>>> xstFreq 250
>>>>> outputEnergies 100
>>>>> outputPressure 100
>>>>>
>>>>>
>>>>> minimize 1000
>>>>> run 50000
>>>>>
>>>>> On Fri, Mar 11, 2016 at 3:34 PM, Maxim Belkin <mbelkin_at_ks.uiuc.edu> wrote:
>>>>>> It is better to set cellbasisvector[1,2] to exact values.
>>>>>> Because you are not using MSM (and you can’t use PME with non-periodic systems), interactions between atoms in the system aren’t right (they are cut off at 12 AA). You should look into MSM first, it should not be that difficult.
>>>>>>
>>>>>> Now, the fact that your simulation crash ~150,000 steps might indicate a problem that emerges due to multiple time-stepping. Try: "stepspercycle 2".
>>>>>>
>>>>>> Finally, what do you expect to get from an NPT simulation with vacuum along z? How do you expect langevinPiston to work in this situation?
>>>>>>
>>>>>>
>>>>>>
>>>>>>> On Mar 10, 2016, at 20:01, Olya Kravchenko <ovkrav_at_gmail.com> wrote:
>>>>>>>
>>>>>>> Thank you! I will study MSM method.
>>>>>>>
>>>>>>> I have changed the size of cell vectors back to 62, that is the actual
>>>>>>> size rounded up, and I commented out cellBasisVector3. I ran a number
>>>>>>> of simulations today and I keep seeing the same error when I increase
>>>>>>> the number of timesteps. The margin size does seem to make the job run
>>>>>>> longer, at least I don't see any other correlation. My last run quit
>>>>>>> near 150 000 timesteps (I set it up for 500K timesteps).
>>>>>>>
>>>>>>> Apart from that, I see in vmd that my simulation box actually becomes
>>>>>>> longer in Z-direction, it is a cube in the beginning and then the
>>>>>>> proportion changes. I am not sure how to interpret it, but it looks
>>>>>>> like the size in z does not hold, certainly no hard wall. I wonder if
>>>>>>> there is a way to set up the hard wall in a simple way?
>>>>>>>
>>>>>>> Olga
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> On Thu, Mar 10, 2016 at 4:58 PM, Maxim Belkin <mbelkin_at_ks.uiuc.edu> wrote:
>>>>>>>> Olya,
>>>>>>>>
>>>>>>>> If your system size is 60A in x and y, you should use vectors of length 60 as your "cellbasisvector[1,2]". You are currently using 70, which means you have vacuum in x and y directions and this is not what you want. To simulate a system that is non-periodic along Z with NAMD you have to use Multi-level summation method (look for MSM option in NAMD User’s Guide). I believe you have to remove "cellbasisvector3" instead of specifying 0 0 0 but I’m not 100% positive on that.
>>>>>>>>
>>>>>>>> Maxim
>>>>>>>>
>>>>>>>>
>>>>>>>>> On Mar 10, 2016, at 11:21, Olya Kravchenko <ovkrav_at_gmail.com> wrote:
>>>>>>>>>
>>>>>>>>> My system consists of 200 argon atoms in a box, PBC in X and Y, I
>>>>>>>>> would like hard wall in Z and my goal is to measure
>>>>>>>>> density/concentration. So I think that NPT is reasonable in my case,
>>>>>>>>> is that correct?
>>>>>>>>>
>>>>>>>>> I commented out PME, because I read that it will speed up
>>>>>>>>> calculations, also because, if I understand it correctly, it requires
>>>>>>>>> all three cell basis vectors to be defined.
>>>>>>>>>
>>>>>>>>> linaccel is commented because I defined my own force in a separate
>>>>>>>>> file for the force.
>>>>>>>>>
>>>>>>>>> Here is my configuration file:
>>>>>>>>>
>>>>>>>>> #############################################################
>>>>>>>>> ## JOB DESCRIPTION ##
>>>>>>>>> #############################################################
>>>>>>>>>
>>>>>>>>> # Equilibration of argon atoms with PBC in XY direction
>>>>>>>>>
>>>>>>>>> #############################################################
>>>>>>>>> ## ADJUSTABLE PARAMETERS ##
>>>>>>>>> #############################################################
>>>>>>>>>
>>>>>>>>> structure ar.psf
>>>>>>>>> coordinates ar.pdb
>>>>>>>>>
>>>>>>>>> set temperature 310
>>>>>>>>> set outputname ar_force_out
>>>>>>>>>
>>>>>>>>> firsttimestep 0
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> #############################################################
>>>>>>>>> ## SIMULATION PARAMETERS ##
>>>>>>>>> #############################################################
>>>>>>>>>
>>>>>>>>> # Input
>>>>>>>>> paraTypeCharmm on
>>>>>>>>> parameters My_parameters.rtf
>>>>>>>>> temperature $temperature
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> # Force-Field Parameters
>>>>>>>>> exclude scaled1-4
>>>>>>>>> 1-4scaling 1.0
>>>>>>>>> cutoff 12.0
>>>>>>>>> switching on
>>>>>>>>> switchdist 10.0
>>>>>>>>> pairlistdist 14.0
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> # Integrator Parameters
>>>>>>>>> timestep 2.0 ;# 2fs/step
>>>>>>>>> rigidBonds all ;# needed for 2fs steps
>>>>>>>>> nonbondedFreq 1
>>>>>>>>> fullElectFrequency 2
>>>>>>>>> stepspercycle 10
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> # Constant Temperature Control
>>>>>>>>> langevin on ;# do langevin dynamics
>>>>>>>>> langevinDamping 1 ;# damping coefficient (gamma) of 1/ps
>>>>>>>>> langevinTemp $temperature
>>>>>>>>> langevinHydrogen off ;# don't couple langevin bath to hydrogens
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> # Periodic Boundary Conditions
>>>>>>>>> cellBasisVector1 62.0 0. 0.0
>>>>>>>>> cellBasisVector2 0.0 62.0 0.0
>>>>>>>>> cellBasisVector3 0.0 0 0.0
>>>>>>>>> cellOrigin 30.0 30.0 32.5
>>>>>>>>>
>>>>>>>>> wrapAll on
>>>>>>>>>
>>>>>>>>> margin 10
>>>>>>>>>
>>>>>>>>> # PME (for full-system periodic electrostatics)
>>>>>>>>> #PME yes
>>>>>>>>> #PMEGridSpacing 1.0
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> # Constant Pressure Control (variable volume)
>>>>>>>>> useGroupPressure yes ;# needed for rigidBonds
>>>>>>>>> useFlexibleCell no
>>>>>>>>> useConstantArea no
>>>>>>>>>
>>>>>>>>> langevinPiston on
>>>>>>>>> langevinPistonTarget 1.01325 ;# in bar -> 1 atm
>>>>>>>>> langevinPistonPeriod 100.0
>>>>>>>>> langevinPistonDecay 50.0
>>>>>>>>> langevinPistonTemp $temperature
>>>>>>>>>
>>>>>>>>> ##################
>>>>>>>>> ###DEFINE FORCES##
>>>>>>>>> ##################
>>>>>>>>>
>>>>>>>>> tclForces on
>>>>>>>>> tclForcesScript my_force.tcl
>>>>>>>>> #set linaccel "30 0 0"
>>>>>>>>>
>>>>>>>>> # Output
>>>>>>>>> outputName $outputname
>>>>>>>>>
>>>>>>>>> restartfreq 500 ;# 500steps = every 1ps
>>>>>>>>> dcdfreq 250
>>>>>>>>> xstFreq 250
>>>>>>>>> outputEnergies 100
>>>>>>>>> outputPressure 100
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> # Minimization
>>>>>>>>> minimize 100
>>>>>>>>> reinitvels $temperature
>>>>>>>>>
>>>>>>>>> run 100000
>>>>>>>>>
>>>>>>>>> Thanks!!
>>>>>>>>>
>>>>>>>>> Olga
>>>>>>>>>
>>>>>>>>> On Thu, Mar 10, 2016 at 2:05 AM, Norman Geist
>>>>>>>>> <norman.geist_at_uni-greifswald.de> wrote:
>>>>>>>>>> We need to see your full input script, since it looks you're doing npt?
>>>>>>>>>>
>>>>>>>>>> Norman Geist
>>>>>>>>>>
>>>>>>>>>>> -----Ursprüngliche Nachricht-----
>>>>>>>>>>> Von: owner-namd-l_at_ks.uiuc.edu [mailto:owner-namd-l_at_ks.uiuc.edu] Im
>>>>>>>>>>> Auftrag von Olya Kravchenko
>>>>>>>>>>> Gesendet: Donnerstag, 10. März 2016 04:23
>>>>>>>>>>> An: NAMD <namd-l_at_ks.uiuc.edu>
>>>>>>>>>>> Betreff: namd-l: periodic boundary condition in XY, hard wall in Z - errors
>>>>>>>>>>>
>>>>>>>>>>> Hi all,
>>>>>>>>>>>
>>>>>>>>>>> I would like to apply periodic boundary conditions only in X and Y
>>>>>>>>>>> directions. In my config file I used the following entry:
>>>>>>>>>>>
>>>>>>>>>>> # Periodic Boundary Conditions
>>>>>>>>>>> cellBasisVector1 70.0 0. 0.0
>>>>>>>>>>> cellBasisVector2 0.0 70.0 0.0
>>>>>>>>>>> cellBasisVector3 0.0 0 0.0
>>>>>>>>>>> cellOrigin 30.0 30.0 32.5
>>>>>>>>>>>
>>>>>>>>>>> wrapAll on
>>>>>>>>>>> margin 5
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> When I submit the job it runs for a while until I see the following error:
>>>>>>>>>>>
>>>>>>>>>>> WRITING COORDINATES TO DCD FILE ar_force_out.dcd AT STEP 20750
>>>>>>>>>>> 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.
>>>>>>>>>>>
>>>>>>>>>>> How can I address this error? Is this the right way to set up PBC in 2
>>>>>>>>>>> dimensions?
>>>>>>>>>>>
>>>>>>>>>>> My simulation cell is actually around (~61 ~61 ~62), as measured by
>>>>>>>>>>> minmax. I played with numbers and noticed that my simulation runs
>>>>>>>>>>> longer if I increase both the margin and cell basis vectors.
>>>>>>>>>>>
>>>>>>>>>>> I saw some of the older threads in the archives where the margin was
>>>>>>>>>>> discussed but it's not clear how large should it be. If I want to run
>>>>>>>>>>> my simulation for 500 000 steps, what margin should I set up?
>>>>>>>>>>>
>>>>>>>>>>> Also, when choosing cell basis vector, how much bigger should it be
>>>>>>>>>>> than what is measured by "measure minmax"?
>>>>>>>>>>
>>>>>>>>
>>>>>>
>>>>
>>

This archive was generated by hypermail 2.1.6 : Sun Dec 31 2017 - 23:20:14 CST