Re: [Towhee-users] VOLNPT output: an error message?

[Aric N. <agn...@en...>]

Dear Marcus,

Thank you for your reply.  I understand that the errant output is a result of the model violating the constraints of the pmcell move by having bonds in each montmorillonite layer that cross the boundary.  The harmonic bonds that cross the boundary have spring constants of zero.  Thus, the energy of these bonds is independent of the positions of the atoms (i.e. it’s always zero).  This condition is a result of the parameters in the CLAYFF forcefield.  I have only “created” these bonds because I would like each montmorillonite layer to be able to translate (as a unit)  in the x-, y-, or z-direction as part of the MC moves.

I can create an input file without these zero-force constant “bonds”, in which case each atom in the montmorillonite is a molecule, but I have been unable to group these individual molecules into a super-molecule that I can then translate as a unit.  My starting simulation cell size is the equilibrium cell size from a 1.0 ns MD-NPT in LAMMPS, so I am not too worried about the changes in the MMT density.  With an anisotropic MC-NPT that only has a five percent probability of performing a pmcell move, I suspect that density fluctuations in the montmorillonite should be limited and, I am assuming, negligible.  Is this a reasonable assumption?

These simulations are part of a thermodynamic integration.  The model has only one box and the equilibrium cell parameters from the MD-NPT simulation confirm that my layer spacing and water content are consistent with the experimental condition that I am trying to simulate.  So, if the pmcell move adjusts the center of mass of each atom in the montmorillonite sheets in a manner appropriate with the volume change (adjusting interatomic distances), then that is what I want.  However, is there a way to also get the center of mass of each sheet to translate as part of the pmtracm command when in fact each sheet is composed of many molecules?  The pmplane move sounds close to what I want, but I am unsure how to define a dynamic definition of the plane so that it always includes only a montmorillonite sheet regardless of previous translation moves during the simulation. Do you have any suggestions or is simulating with rigid montmorillonite layers the only option with the current move algorithms?

Thank you,

Aric Newton


On Oct 6, 2014, at 10:02 AM, Marcus Martin <me...@gm...> wrote:

> On Fri, Oct 3, 2014 at 4:22 AM, Aric Newton <agn...@en...> wrote:
> 
> I am using towhee for a one box MC-NPT simulation of two montmorillonite layers with an interlayer space containing the counter ions and water molecules.  Each montmorillonite layer is treated as an individual molecule and moves as a unit via the pmtracm move probability.  I have previously equilibrated the system at high temperature in the NVT ensemble without any errors, but the output that I am getting now that I have switched to the NPT ensemble leaves me a bit uncertain.  In the towhee_output file, there are about 10,000 lines that read as below.   
> 
>  VOLNPT: move attemped in box that has a periodic bonded molecule
>  Molecule type:      -31414
>  Chain number:           1
>  Box:           2
> ...
> 
> The *.pdb snapshots are being output during the simulation, but the rest of the output appears to be non-standard/ non-sensical- there is only one box in the simulation and the negative number as a molecule type seems odd.
> 
> 
> Indeed the error messages do mean there is a problem.  When you perform a volume move in Towhee it scales the center of mass position of all the molecules (in a manner appropriate for the volume change) and then keeps all of the intramolecular distances the same.  This is not possible when the "molecule" spans the periodic images with its bonding pattern.
> 
> You would need a completely different move that scaled the intramolecular distances for the volume change.  That would then lead to problems with any molecule that has a rigid bond length (as the attempt to change that bond length would result in an infinite energy) and also would be fairly inefficient as the intramolecular distances are quite sensitive.
>  
> The volnpt.f file description from the towhee code manual notes that “you cannot use this move with a molecule that is bonded through the periodic boundary condition (like a zeolite).”  I suspect that the bonds in the montmorillonite molecules are across the periodic boundary and trigger the VOLNPT output message. However, this same note is listed for the volnvt.f file and there are no warnings generated.  Are these simulation results corrupt?  Is there a way to salvage this simulation or work around this problem?  The forcefield for the interactions is CLAYFF.  So, the bonds, angles, torsions, and improper  in the montmorillonite layers are zero energy.  I had to define these terms so that each montmorillonite layer can be a molecule that can move/ translate during the simulation.  Is there an alternate way to conduct this simulation and avoid the VOLNPT violation?
> 
> Only if every atom in the montmorillonite was a separate molecule.  However, be aware you are then changing the density of your montmorillonite and it is quite likely that is not what you are intending to sample.
> 
> I'm not sure if you were intending to sample the layer separation distance, but in that case you would need a new implementation of the volume move that only scale the volume in one of the dimensions instead of uniformly.
> 
> Marcus
> 
> -- 
> Marcus G. Martin
> 88 Martinez Road
> Edgewood NM 87015-8222
> land (505) 286 4457
> cell (505) 363 3179
> www.photobirder.com
> towhee.sourceforge.net
> 

Aric Newton

Laboratory of Nuclear and Environmental Materials
Faculty of Engineering
Hokkaido University
Sapporo, Japan