crystal Modeling

Contributed by Guido Raos for pDynamo version 1.7.2.

The code on this page is available here.

Introduction

pDynamo has several features allowing the simulation of organic crystal structures. These include:

• The possibility of simulating systems with general periodic boundary conditions, using unit cells belonging to different crystal classes (cubic, orthorhombic, monoclinic, triclinic, etc.).

• The possibility of specifying any crystal symmetry through the fractional coordinates of the symmetry-equivalent points within a unit cell. These symmetry information can be used in the construction of the initial crystal structure, and if desired it can be maintained during the subsequent simulation steps (i.e., optimization of the crystal structure).

• The possibility of simulating "small" crystal cells, with lengths smaller than twice the cutoff for the non-bonded interactions (typically, 2.5-3.0 nm). This is possible because pDynamo does not use the common "minimum image" convention.

Some of these features will be illustrated below using paracetamol as an example.

Example: modeling paracetamol crystals

The definition of a force field for paracetamol is a prerequisite for this tutorial, which is covered in another section of this site. Below, we will assume that a ParacetamolMM-0.yaml file and a ParacetamolQM-0.yaml file have been produced and tested by standard gas-phase calculations, as described in that section (the QM method used for producing the second yaml file is PM6).

The starting structure for our simulations will be the HXACAN19 entry in the Cambridge Structural Database (CSD). Typically, this would be available as a Crystallographic Information File (HXACAN19.cif). pDynamo does not have the capability of reading information directly from such files, yet. At the moment, the user has to convert the cif file into another pDynamo-friendly file (mol or pdb, for example), process this file to read in the atomic coordinates and define the parameters for the MM or QM calculations (these steps will be skipped here by using the Paracetamol*-0.yaml files), and then re-insert "manually" the crystallographic information, as explained below. Notice also that, at the moment, pDynamo's PDBFile_ToSystem() and PDBFile_FromSystem() do not read/write any crystallographic information in the CRYST1 record of PDB files.

The most important crystallographic information within the HXACAN19.cif file is that paracetamol crystallizes within the monoclinic crystal system, in space group P2₁/a. The lattice parameters characterizing the unit cell are: a=12.872 Å, b=9.370 Å, c=7.085 Å, β=115.62° (this being a monoclinic unit cell, it is implicit that α=γ=90°).

The symmetry elements of a P2₁/a cell include a screw axis, a glide plane and an inversion center, as well as their combinations. This symmetry information can be expressed as follows: for each point (in particular, this might be the location of an atom) within the unit cell with fractional coordinates (x,y,z), there are three other points related to it by symmetry with fractional coordinates (1/2-x,1/2+y,-z), (-x,-y,-z) and (1/2+x,1/2-y,z). This type of symmetry information can be encoded within a Python dictionary, which will then be imported by the pDynamo script doing the actual simulations (it should be very easy to extend this list of space groups, whenever there is need of a new type of crystal symmetry).

SpaceGroupDefinitions.py:

"""A handy collection of crystal-related definitions for pDynamo"""

from pCore import Transformation3_FromSymmetryOperationString, Transformation3Container

from pMolecule import CrystalClassCubic, CrystalClassHexagonal, CrystalClassMonoclinic, CrystalClassOrthorhombic, \

CrystalClassRhombohedral, CrystalClassTetragonal, CrystalClassTriclinic

from pMoleculeScripts import CrystalAnalyzeTransformations, CrystalExpandToP1, \

CrystalGetImageBondPairs, CrystalCenterCoordinates

cubic = CrystalClassCubic()

hexagonal = CrystalClassHexagonal()

monoclinic = CrystalClassMonoclinic()

orthorhombic = CrystalClassOrthorhombic()

rhombohedral = CrystalClassRhombohedral()

tetragonal = CrystalClassTetragonal()

triclinic = CrystalClassTriclinic()

# . Python dictionary containing the crystal class and the symmetry operations for a few space groups

spacegroup = {

"P1" : [ triclinic, "(x,y,z)" ],

"P-1" : [ triclinic, "(x,y,z)", "(-x,-y,-z)" ],

"C2" : [ monoclinic, "(x,y,z)", "(-x,y,-z)"],

"P21c" : [ monoclinic, "(x,y,z)", "(-x,y+1/2,-z+1/2)", "(-x,-y,-z)", "(x,-y+1/2,z+1/2)" ],

"P21a" : [ monoclinic, "(x,y,z)", "(-x+1/2,y+1/2,-z)", "(-x,-y,-z)", "(x+1/2,-y+1/2,z)" ],

"P 1 21/A 1" : [ monoclinic, "(x,y,z)", "(-x+1/2,y+1/2,-z)", "(-x,-y,-z)", "(x+1/2,-y+1/2,z)" ],

}

The pictures below illustrate the symmetry operations of the P2₁/a space group, by showing their action on the unit cell of paracetamol (views along the b and c axes, hydrogens are omitted for clarity):

The following script illustrates the capabilities of pDynamo (see point 1)-3) above), by carrying out the energy evaluation on an MM model of the paracetamol unit cell, both with and without symmetry constraint, as well as a supercell constructed by replicating the original unit cell a few times along x, y and z:

CrystalMMEnergy.py:

"""Definition of crystal symmetry and simple energy calculations on Paracetamol by MM"""

from pCore import logFile, Clone, YAMLPickle, YAMLUnpickle

from pMolecule import NBModelGABFS

from SpaceGroupDefinitions import *

# . Read in the system definition and re-calculate the energy from the gas-phase optimization.

# . Redefine the NB model from Full to GABFS (we truncate the non-bonded interactions).

molecule = YAMLUnpickle( "ParacetamolMM-0.yaml" )

nbModel = NBModelGABFS( )

molecule.DefineNBModel( nbModel )

gasEnergy = molecule.Energy()

# . Define the crystal by giving the space group information and the cell parameters

crystal = Clone( molecule )

crystal.label = "Paracetamol MM crystal, unit cell with P21/a symmetry"

name = "P21a"

crystalclass = spacegroup[name][0]

trans = []

for t in spacegroup[name][1:]:

trans.append( Transformation3_FromSymmetryOperationString(t) )

transcontainer = Transformation3Container.WithTransformations( trans )

crystal.DefineSymmetry( crystalClass = crystalclass, transformations=transcontainer, \

a=12.872, b=9.370, c=7.085, beta=115.62)

crystal.symmetry.Summary( )

CrystalAnalyzeTransformations( crystal )

zfactor = float( len(crystal.symmetry.transformations) )

# . Obtain some information about the system

cellParameter = crystal.symmetry.crystalClass.GetUniqueSymmetryParameters( crystal.symmetryParameters )

# . Calculate the energy and save to file for later use

energy = crystal.Energy() * zfactor

YAMLPickle( "ParacetamolMM-P21a.yaml", crystal )

# . Define a new system in which (almost) all symmetry is removed

p1crystal = CrystalExpandToP1( crystal )

p1crystal.label = "Paracetamol MM crystal, unit cell with P1 symmetry"

p1cellParameter = \

p1crystal.symmetry.crystalClass.GetUniqueSymmetryParameters( p1crystal.symmetryParameters )

p1energy = p1crystal.Energy()

YAMLPickle( "ParacetamolMM-P1.yaml", p1crystal )

# . Define a new system made up by a big supercell

na=3; nb=3; nc=5; ncells=na*nb*nc

bigcrystal = Clone( p1crystal )

bigcrystal.label = "Paracetamol MM crystal, big supercell"

bigcrystal = CrystalExpandToP1( bigcrystal, arange=range(na), brange=range(nb), crange=range(nc) )

bigcellParameter = \

bigcrystal.symmetry.crystalClass.GetUniqueSymmetryParameters( bigcrystal.symmetryParameters )

bigenergy = bigcrystal.Energy() / ncells

# . Print the results.

table = logFile.GetTable( columns = [ 22, 10, 10, 10, 10, 22, 22 ] )

table.Start( )

table.Title( "Comparison of Crystal Energies" )

table.Heading( "System description" )

table.Heading( "a" ); table.Heading ( "b" ); table.Heading ( "c" ); table.Heading ( "beta" );

table.Heading( "Tot.Energy/cell" )

table.Heading( "Subl.Energy/mol" )

table.Entry( "Unit cell (P21/a)", alignment = "l" )

table.Entry( "%8.4f" % cellParameter['a'] )

table.Entry( "%8.4f" % cellParameter['b'] )

table.Entry( "%8.4f" % cellParameter['c'] )

table.Entry( "%8.4f" % cellParameter['beta'] )

table.Entry( "%20.8f" % ( energy ) )

table.Entry( "%20.8f" % ( gasEnergy - energy/zfactor ) )

table.Entry( "Unit cell (P1)", alignment = "l" )

table.Entry( "%8.4f" % p1cellParameter['a'] )

table.Entry( "%8.4f" % p1cellParameter['b'] )

table.Entry( "%8.4f" % p1cellParameter['c'] )

table.Entry( "%8.4f" % p1cellParameter['beta'] )

table.Entry( "%20.8f" % ( p1energy ) )

table.Entry( "%20.8f" % ( gasEnergy-p1energy/zfactor ) )

cellDescription = "Supercell-" + str(na) + "a" + str(nb) + "b" + str(nc) + "c"

table.Entry( cellDescription, alignment = "l" )

table.Entry( "%8.4f" % bigcellParameter['a'] )

table.Entry( "%8.4f" % bigcellParameter['b'] )

table.Entry( "%8.4f" % bigcellParameter['c'] )

table.Entry( "%8.4f" % bigcellParameter['beta'] )

table.Entry( "%20.8f" % bigenergy )

table.Entry( "%20.8f" % ( gasEnergy-bigenergy/zfactor ) )

table.Stop( )

The script illustrates several functionalities of pDynamo, including:

• • the use of the DefineSymmetry method, which expects as input the crystal class, the space group (defined in terms of a "transformation container"), and the lattice parameters;

  • The CrystalAnalyzeTransformations function, which analyses the symmetry transformation in terms of combinations of translations and rotations (which in turn can be proper or improper, the latter corresponding to a reflection through a mirror plane or an inversion center);

  • the use of the GetUniqueSymmetryParameters method to re-extract information about the lattice parameters (this is really unnecessary now, but it will be useful later, once we start re-optimizing the cell);

  • the CrystalExpandToP1 function, which can remove the symmetry constraints on a unit cell (P1 is the space group corresponding to "no symmetry, except for translation"), but it can also create a larger supercell when the na, nb and nc arguments are larger than unity.

The final table summarizing the results from CrystalMMEnergy.py is:

...

----------------------------------------------------------------------------------------------------------

Comparison of Crystal Energies

----------------------------------------------------------------------------------------------------------

System description a b c beta Tot.Energy/cell Subl.Energy/mol

----------------------------------------------------------------------------------------------------------

Unit cell (P21/a) 12.8720 9.3700 7.0850 115.6200 -279.77155841 136.18401375

Unit cell (P1) 12.8720 9.3700 7.0850 115.6200 -279.77155841 136.18401375

Supercell-3a3b5c 38.6160 28.1100 35.4250 115.6200 -279.77155841 136.18401375

----------------------------------------------------------------------------------------------------------

Notice that, to make the energies comparable, we had to multiply/divide them by an appropriate factor. Having done so, all the energies are exactly identical, as it should be since the energy is an extensive quantity and we are working with an additive MM model.

An analogous calculation on the same crystal structure, but using the semiempirical PM6 Hamiltonian, would be:

CrystalQMEnergy.py:

"""Definition of crystal symmetry and simple energy calculations on Paracetamol by QM"""

from pCore import logFile, Clone, YAMLPickle, YAMLUnpickle

from pMolecule import NBModelGABFS, QCModelMNDO

from SpaceGroupDefinitions import *

nbModel = NBModelGABFS( )

qcModel = QCModelMNDO( "pm6" )

# . Read in the system definition and re-calculate the energy from the gas-phase optimization.

# . Redefine the NB model from Full to GABFS (we truncate the non-bonded interactions).

molecule = YAMLUnpickle( "ParacetamolQM-0.yaml" )

molecule.DefineNBModel( nbModel )

gasEnergy = molecule.Energy()

# . Define the crystal by giving the space group information and the cell parameters

crystal = Clone( molecule )

crystal.label = "Paracetamol QM crystal, unit cell with P21/a symmetry"

name = "P21a"

crystalclass = spacegroup[name][0]

trans = []

for t in spacegroup[name][1:]:

trans.append( Transformation3_FromSymmetryOperationString(t) )

transcontainer = Transformation3Container.WithTransformations( trans )

crystal.DefineSymmetry( crystalClass = crystalclass, transformations=transcontainer, \

a=12.872, b=9.370, c=7.085, beta=115.62)

crystal.symmetry.Summary( )

CrystalAnalyzeTransformations( crystal )

zfactor = float( len(crystal.symmetry.transformations) )

# . Obtain some information about the system

cellParameter = crystal.symmetry.crystalClass.GetUniqueSymmetryParameters( crystal.symmetryParameters )

# . Calculate the energy and save to file for later use

crystal.DefineQCModel( qcModel )

crystal.DefineNBModel( nbModel )

energy = crystal.Energy() * zfactor

YAMLPickle( "ParacetamolQM-P21a.yaml", crystal )

# . Define a new system in which (almost) all symmetry is removed

p1crystal = CrystalExpandToP1( crystal )

p1crystal.label = "Paracetamol QM crystal, unit cell with P1 symmetry"

p1cellParameter = \

p1crystal.symmetry.crystalClass.GetUniqueSymmetryParameters( p1crystal.symmetryParameters )

p1crystal.DefineQCModel( qcModel )

p1crystal.DefineNBModel( nbModel )

p1energy = p1crystal.Energy()

YAMLPickle( "ParacetamolQM-P1.yaml", p1crystal )

# . Define a new system made up by a not-so-big supercell

na=1; nb=2; nc=2; ncells=na*nb*nc

bigcrystal = Clone( p1crystal )

bigcrystal.label = "Paracetamol QM crystal, big supercell"

bigcrystal = CrystalExpandToP1( bigcrystal, arange=range(na), brange=range(nb), crange=range(nc) )

bigcellParameter = \

bigcrystal.symmetry.crystalClass.GetUniqueSymmetryParameters( bigcrystal.symmetryParameters )

bigcrystal.DefineQCModel( qcModel )

bigcrystal.DefineNBModel( nbModel )

bigenergy = bigcrystal.Energy() / ncells

# . Print the results.

table = logFile.GetTable( columns = [ 22, 10, 10, 10, 10, 22, 22 ] )

table.Start( )

table.Title( "Comparison of Crystal Energies" )

table.Heading( "System description" )

table.Heading( "a" ); table.Heading ( "b" ); table.Heading ( "c" ); table.Heading ( "beta" );

table.Heading( "Tot.Energy/cell" )

table.Heading( "Subl.Energy/mol" )

table.Entry( "Unit cell (P21/a)", alignment = "l" )

table.Entry( "%8.4f" % cellParameter['a'] )

table.Entry( "%8.4f" % cellParameter['b'] )

table.Entry( "%8.4f" % cellParameter['c'] )

table.Entry( "%8.4f" % cellParameter['beta'] )

table.Entry( "%20.8f" % ( energy ) )

table.Entry( "%20.8f" % ( gasEnergy - energy/zfactor ) )

table.Entry( "Unit cell (P1)", alignment = "l" )

table.Entry( "%8.4f" % p1cellParameter['a'] )

table.Entry( "%8.4f" % p1cellParameter['b'] )

table.Entry( "%8.4f" % p1cellParameter['c'] )

table.Entry( "%8.4f" % p1cellParameter['beta'] )

table.Entry( "%20.8f" % ( p1energy ) )

table.Entry( "%20.8f" % ( gasEnergy-p1energy/zfactor ) )

cellDescription = "Supercell-" + str(na) + "a" + str(nb) + "b" + str(nc) + "c"

table.Entry( cellDescription, alignment = "l" )

table.Entry( "%8.4f" % bigcellParameter['a'] )

table.Entry( "%8.4f" % bigcellParameter['b'] )

table.Entry( "%8.4f" % bigcellParameter['c'] )

table.Entry( "%8.4f" % bigcellParameter['beta'] )

table.Entry( "%20.8f" % bigenergy )

table.Entry( "%20.8f" % ( gasEnergy-bigenergy/zfactor ) )

table.Stop( )

Notice that all the atoms within the asymmetric unit (this coincides with one molecule under P2₁/a, one unit cell under P1) are treated at the QM level. Instead, the interaction between this asymmetric unit and the remainder of the crystal is treated at the MM level (i.e., as a sum of Lennard-Jones and charge-charge interactions, since there are no link atoms in this system). Notice also that, in the last part of the calculation, we have used a small 1x2x2 supercell to reduce computer time.

The last part of the output from CrystalQMEnergy.py is now:

...

----------------------------------------------------------------------------------------------------------

Comparison of Crystal Energies

----------------------------------------------------------------------------------------------------------

System description a b c beta Tot.Energy/cell Subl.Energy/mol

----------------------------------------------------------------------------------------------------------

Unit cell (P21/a) 12.8720 9.3700 7.0850 115.6200 -1206.32319623 121.72420222

Unit cell (P1) 12.8720 9.3700 7.0850 115.6200 -1173.20622249 113.44495878

Supercell-1a2b2c 12.8720 18.7400 14.1700 115.6200 -1115.30385304 98.96936642

----------------------------------------------------------------------------------------------------------

Notice that, unlike the MM case, the energies are not independent of system size. This is not due to an error in the program. Rather, it depends on two factors: (a) QM models include polarization (i.e., many body) effects and therefore they are not rigorously additive, and (b) as we expand the asymmetric unit, we change the relative weights of MM-MM, QM-MM and QM-QM non-bonded interactions.

Finally, the following script illustrates a cascade of MM minimizations of the paracetamol crystal structure. We start from the original P21/a unit cell and progressively remove the constraints on it (these include both restrictions on the symmetry of the system and on the values of the lattice parameters). The aim is essentially to check the energetic stability of the crystal structure and/or the quality of the force field employed in the simulation.

CrystalMMMinimizations.py:

"""Example 2 for crystal structures: cascade of energy minimizations on Paracetamol by MM"""

# . Script by guido.raos@polimi.it, tested with pDynamo-1.7.2 (September 2012)

from Definitions import *

from pCore import YAMLPickle, YAMLUnpickle

from pBabel import PDBFile_FromSystem

from SpaceGroupDefinitions import *

# . Initialize some useful lists

crystal = []

description = []

cellParam = []

energy = []

rmsd = []

# . Starting system

i = 0

description.append( "Unit cell 0 (P21/a)" )

logFile.Paragraph( "CRYSTAL " + str(i) + ": " + description[i] )

crystal.append( YAMLUnpickle( "ParacetamolMM-P21a.yaml" ) )

cp = crystal[i].symmetry.crystalClass.GetUniqueSymmetryParameters ( crystal[i].symmetryParameters )

cp['alpha'] = 90.0; cp['gamma'] = 90.0

cellParam.append( cp )

zfactor = float( len(crystal[0].symmetry.transformations) )

energy.append( crystal[i].Energy() * zfactor )

rmsd.append( 0.0 )

# . Define a few useful quantities before proceeding to the other systems

#(coord0, translation) = CrystalCenterCoordinates( crystal[i] )

coord0 = Clone( crystal[0].coordinates3 )

masses = crystal[0].atoms.GetItemAttributes( "mass" )

crystal0p1 = CrystalExpandToP1( crystal[i] )

#(coord0p1, translation) = CrystalCenterCoordinates( crystal0p1 )

coord0p1 = Clone( crystal0p1.coordinates3 )

massesp1 = crystal0p1.atoms.GetItemAttributes( "mass" )

# . Optimization with fixed lattice parameters and full P21/a symmetry

i += 1

description.append( "Unit cell 1 (P21/a)" )

logFile.Paragraph( "CRYSTAL " + str(i) + ": " + description[i] )

crystal.append( Clone( crystal[i-1] ) )

ConjugateGradientMinimize_SystemGeometry ( crystal[i], logFrequency=10, maximumIterations=1000, \

rmsGradientTolerance=0.5 )

cp = crystal[i].symmetry.crystalClass.GetUniqueSymmetryParameters ( crystal[i].symmetryParameters )

cp['alpha'] = 90.0; cp['gamma'] = 90.0

cellParam.append( cp )

energy.append( crystal[i].Energy() * zfactor )

#(coord3, translation) = CrystalCenterCoordinates( crystal[i] )

coord3 = Clone( crystal[i].coordinates3 )

coord3.Superimpose( coord0, weights = masses )

rmsd.append( coord3.RMSDeviation( coord0, weights = masses ) )

# . Optimization with variable lattice parameters and full P21/a symmetry

i += 1

description.append( "Unit cell 2 (P21/a)" )

logFile.Paragraph( "CRYSTAL " + str(i) + ": " + description[i] )

crystal.append( Clone( crystal[i-1] ) )

ConjugateGradientMinimize_SystemGeometry( crystal[i], logFrequency=10, maximumIterations=1000, \

rmsGradientTolerance=0.5, QOPTIMIZESYMMETRYPARAMETERS=True )

cp = crystal[i].symmetry.crystalClass.GetUniqueSymmetryParameters ( crystal[i].symmetryParameters )

cp['alpha'] = 90.0; cp['gamma'] = 90.0

cellParam.append( cp )

energy.append( crystal[i].Energy() * zfactor )

#(coord3, translation) = CrystalCenterCoordinates( crystal[i] )

coord3 = Clone( crystal[i].coordinates3 )

coord3.Superimpose( coord0, weights = masses )

rmsd.append( coord3.RMSDeviation( coord0, weights = masses ) )

# . Optimization with variable lattice parameters and P1 symmetry, but retaining the monoclinic crystal class

i += 1

description.append( "Unit cell 3 (P1, monocl.)" )

logFile.Paragraph( "CRYSTAL " + str(i) + ": " + description[i] )

crystal.append( CrystalExpandToP1( crystal[i-1] ) )

ConjugateGradientMinimize_SystemGeometry( crystal[i], logFrequency=10, maximumIterations=1000, \

rmsGradientTolerance=0.5, QOPTIMIZESYMMETRYPARAMETERS=True )

cp = crystal[i].symmetry.crystalClass.GetUniqueSymmetryParameters ( crystal[i].symmetryParameters )

cp['alpha'] = 90.0; cp['gamma'] = 90.0

cellParam.append( cp )

energy.append( crystal[i].Energy() )

#(coord3, translation) = CrystalCenterCoordinates( crystal[i] )

coord3 = Clone( crystal[i].coordinates3 )

coord3.Superimpose( coord0p1, weights = massesp1 )

rmsd.append( coord3.RMSDeviation( coord0p1, weights = massesp1 ) )

# . Optimization with variable lattice parameters and P1 symmetry, but switching to the triclinic crystal class.

# . QTRICLINIC option of CrystalExpandToP1() has not yet been implemented. We'll do this "manually".

i += 1

description.append( "Unit cell 4 (P1, tricl.)" )

logFile.Paragraph( "CRYSTAL " + str(i) + ": " + description[i] )

crystal.append( Clone( crystal[i-1] ) )

# . Note the slightly "off-90 degrees" values of alpha and gamma.

crystal[i].DefineSymmetry( crystalClass=CrystalClassTriclinic(), \

a=cellParam[i-1]['a'], b=cellParam[i-1]['b'], c=cellParam[i-1]['c'], \

alpha=90.01, beta=cellParam[i-1]['beta'], gamma=89.99 )

ConjugateGradientMinimize_SystemGeometry( crystal[i], logFrequency=10, maximumIterations=1000, \

rmsGradientTolerance=0.5, QOPTIMIZESYMMETRYPARAMETERS=True )

cellParam.append( crystal[i].symmetry.crystalClass.GetUniqueSymmetryParameters( crystal[i].symmetryParameters ) )

energy.append( crystal[i].Energy() )

#(coord3, translation) = CrystalCenterCoordinates( crystal[i] )

coord3 = Clone( crystal[i].coordinates3 )

coord3.Superimpose( coord0p1, weights = massesp1 )

rmsd.append( coord3.RMSDeviation( coord0p1, weights = massesp1 ) )

# . Optimization with lower translational symmetry, by expanding to a bigger P1 supercell.

na=3; nb=3; nc=5; ncells=na*nb*nc

i += 1

description.append( "Supercell-" + str(na) + "a" + str(nb) + "b" + str(nc) + "c" )

logFile.Paragraph( "CRYSTAL " + str(i) + ": " + description[i] )

# . CrystalExpandToP1() does not return a new instance of the system, if this is already P1.

crystal.append( Clone( crystal[i-1] ) )

crystal[i] = CrystalExpandToP1( crystal[i], arange=range(na), brange=range(nb), crange=range(nc))

ConjugateGradientMinimize_SystemGeometry( crystal[i], logFrequency=10, maximumIterations=1000, \

rmsGradientTolerance=0.5, QOPTIMIZESYMMETRYPARAMETERS=True )

cellParam.append( crystal[i].symmetry.crystalClass.GetUniqueSymmetryParameters( crystal[i].symmetryParameters ) )

cellParam[i]["a"] /= na;

cellParam[i]["b"] /= nb;

cellParam[i]["c"] /= nc;

energy.append( crystal[i].Energy()/ncells )

rmsd.append (-1.0 )

# . Summarize the results in a table.

table = logFile.GetTable ( columns = [ 30, 10, 10, 10, 10, 10, 10, 10, 22 ] )

table.Start ( )

table.Title ( "Comparison of Crystal Energy Minimizations for Paracetamol" )

table.Heading ( "System description" )

table.Heading ( "a/na" ); table.Heading ( "b/nb" ); table.Heading ( "c/nc" );

table.Heading( "alpha" ); table.Heading ( "beta" ); table.Heading( "gamma" );

table.Heading ( "RMSD");

table.Heading ( "Tot.Energy/cell" )

for i in range(len(description)):

table.Entry ( description[i], alignment = "l" )

table.Entry ( "%8.4f" % cellParam[i]['a'] )

table.Entry ( "%8.4f" % cellParam[i]['b'] )

table.Entry ( "%8.4f" % cellParam[i]['c'] )

table.Entry ( "%8.4f" % cellParam[i]['alpha'] )

table.Entry ( "%8.4f" % cellParam[i]['beta'] )

table.Entry ( "%8.4f" % cellParam[i]['gamma'] )

table.Entry ( "%8.4f" % rmsd[i] )

table.Entry ( "%20.8f" % energy[i] )

table.Stop ( )

# . Write every structure in PDB format for visualization (we should add the CRYST1 field to PDBFileWriter!)

for i in range(len(description)):

PDBFile_FromSystem( "Crystal"+str(i)+".pdb", crystal[i] )

# . This supercell could be used to start an MD simulation of the crystal

# . (currently only NVE/NVT with fixed lattice params, using VelocityVerletDynamics_SystemGeometry() ).

YAMLPickle( "Paracetamol"+description[5]+".yaml", crystal[5] )

The final part of the output from CrystalMMMinimizations.py is:

...

--------------------------------------------------------------------------------------------------------------------------

Comparison of Crystal Energy Minimizations for Paracetamol

--------------------------------------------------------------------------------------------------------------------------

System description a/na b/nb c/nc alpha beta gamma RMSD Tot.Energy/cell

--------------------------------------------------------------------------------------------------------------------------

Unit cell 0 (P21/a) 12.8720 9.3700 7.0850 90.0000 115.6200 90.0000 0.0000 -279.77155841

Unit cell 1 (P21/a) 12.8720 9.3700 7.0850 90.0000 115.6200 90.0000 0.0727 -637.48628209

Unit cell 2 (P21/a) 12.6341 9.1756 7.1089 90.0000 115.6640 90.0000 0.0730 -642.29395914

Unit cell 3 (P1, monocl.) 12.6341 9.1756 7.1089 90.0000 115.6640 90.0000 0.1306 -642.29395914

Unit cell 4 (P1, tricl.) 12.4843 9.6852 7.0482 89.9891 118.7800 89.9953 0.2611 -644.63721387

Supercell-3a3b5c 12.5023 9.6829 7.0504 90.0026 118.9255 90.0057 -1.0000 -644.65911991

--------------------------------------------------------------------------------------------------------------------------

Notice that the largest change in energy occurs during the first minimization (P2₁/a cell, fixed lattice parameters). Most of this is probably associated with small adjustments in the positions of the H atoms, which are difficult to locate precisely by X-ray diffraction. Subsequent minimizations produce very small changes in the energy, lattice parameters and atomic coordinates, confirming the stability of the experimental crystal structure and the good quality of the force field employed in the simulations.

Missing features

There are a number of crystal-related features which are still missing from the current distribution of pDynamo (1.7.2). These include:

• • Capability of reading/writing crystallographic information from/to CIF or PDB files.

  • Treatment of periodic bonds across unit cells. This would be useful for modeling synthetic polymer crystals (e.g., stereoregular polyolefins, polyamides, etc.), as well as fibrous biopolymers such as some proteins (e.g. collagen) and DNA.

  • NPT MD simulations (variable lattice parameters, fully anisotropic) with non-orthogonal unit cells. Currently, MD simulations on non-orthogonal cells are restricted to the NVT ensemble (fixed lattice parameters).

  • “Exact” methods for the long-range electrostatic interactions (Ewald, etc.) with general (non-orthogonal) periodic boundary conditions.

There are plans to introduce these capabilities in future versions of the program.