Filling With Packmol
While the MDMC universe.fill method is serviceable for simpler systems, it isn’t particularly sophisticated. For more complex filled universes, MDMC has an integration with packmol, a package for creating atomic configurations such as bilayers or mixtures. If you are not aware of how packmol works, then please read the user guide here.
Note to run this tutorial, you need to have Packmol installed; this is already installed in the Docker container. The MDMC interface to packmol is contained in the objects MDMC.MD.packmol.PackmolSetup and MDMC.MD.packmol.PackmolFiller.
[1]:
from MDMC.MD.packmol import PackmolSetup, PackmolFiller
Creating mixtures
The first example shown is a simple system of water mixed in with ethanol. In a 1:1 molecular ratio.
Firstly, we create the molecules needed for the system (See the notebook on creating atomic configurations for information on this.)
(NB: It is highly recommend to add any atom labels, dihedrals, bond angles, etc. to the Molecules at this stage to ensure they are present in the final universe, rather than having to go in and add them all manually later.)
[2]:
from MDMC.MD import Atom, Molecule, Bond, BondAngle
H1 = Atom('H')
H2 = Atom('H', position=[0., 1.63298, 0.])
O = Atom('O', position=[0., 0.81649, 0.57736])
h2o_bonds = Bond((H1, O), (H2, O))
HOH_angle = BondAngle(H1, O, H2)
water_molecule = Molecule(atoms=[H1, H2, O], interactions=[HOH_angle, h2o_bonds], name="water")
h_1 = Atom("H", position=[-1.9237, 0.3850, 0.0000])
h_2 = Atom("H", position=[2.0985, 0.2306, 0.0000])
h_3 = Atom("H", position=[1.1184, -1.0093, 0.8869])
h_4 = Atom("H", position=[1.1184, -1.0093, -0.8869])
h_5 = Atom("H", position=[-0.0227, 1.1812, 0.8852])
h_6 = Atom("H", position=[-0.0227, 1.1812, -0.8852])
c_1 = Atom("C", position=[1.1879, -0.3829, 0.0000])
c_2 = Atom("C", position=[0.0000, 0.5526, 0.0000])
o_1 = Atom("O", position=[-1.1867, -0.2472, 0.0000])
ch_bond = Bond((h_2, c_1), (h_3, c_1), (h_4, c_1), (h_5, c_2), (h_6, c_2))
co_bond = Bond((c_2, o_1))
cc_bond = Bond((c_1, c_2))
oh_bond = Bond((o_1, h_1))
oh_bond_angle = BondAngle(h_1, o_1, c_2)
ethanol_molecule = Molecule(atoms=[h_1,h_2,h_3,h_4,h_5,h_6,c_1,c_2,o_1], interactions=[ch_bond, co_bond, cc_bond, oh_bond, oh_bond_angle], name="ethanol")
Next, put this into a PackmolSetup object.
Currently, only single fixed molecules, or boxes, cubes and spheres of molecules are supported.
Here we add a cube of molecules with the setup.add_cube method. We provide a desired density in \(Ang^-3\); alternatively one could provide a parameter n_structures for the desired number of Molecules one would like in the space. Later we will see setup.add_box for a non-cubic box; furthermore we have setup.add_sphere to add a sphere of given radius.
[3]:
setup = PackmolSetup()
# Two identically overlapping cubes will create a mixture of water and ethanol in a 1:1 ratio of molecules
setup.add_cube(structure=water_molecule, size=40., density=0.01)
setup.add_cube(structure=ethanol_molecule, size=40., density=0.01)
Next step is to pass this PackmolSetup object to a PackmolFiller object.
Calling fill_with_packmol will setup & run packmol in the background and return a filled universe with the molecules specified.
[4]:
filler = PackmolFiller(setup_data=setup)
universe = filler.fill_with_packmol()
################################################################################
PACKMOL - Packing optimization for the automated generation of
starting configurations for molecular dynamics simulations.
Version 20.14.2
################################################################################
Packmol must be run with: packmol < inputfile.inp
Userguide at: http://m3g.iqm.unicamp.br/packmol
Reading input file... (Control-C aborts)
Types of coordinate files specified: pdb
Seed for random number generator: 1234567
Output file: /home/runner/work/MDMCv0.2_pilot/MDMCv0.2_pilot/doc/how-to/use-MDMC/notebooks/packmol_files/output-universe.pdb
Reading coordinate file: water-0.pdb
Reading coordinate file: ethanol-1.pdb
Number of independent structures: 2
The structures are:
Structure 1 :water-0.pdb( 3 atoms)
Structure 2 :ethanol-1.pdb( 9 atoms)
Maximum number of GENCAN loops for all molecule packing: 400
Total number of restrictions: 2
Distance tolerance: 2.0000000000000000
Warning: Type of residue numbering not set for structure 1
Residue numbering set for structure 1 : 0
Swap chains of molecules of structure 1 : F
Warning: Type of residue numbering not set for structure 2
Residue numbering set for structure 2 : 0
Swap chains of molecules of structure 2 : F
Number of molecules of type 1 : 640
Number of molecules of type 2 : 640
Total number of atoms: 7680
Total number of molecules: 1280
Number of fixed molecules: 0
Number of free molecules: 1280
Number of variables: 7680
Total number of fixed atoms: 0
Maximum internal distance of type 1 : 1.6330000000000000
Maximum internal distance of type 2 : 4.0259464725701459
All atoms must be within these coordinates:
x: [ -1000.0000000000000 , 1000.0000000000000 ]
y: [ -999.20000000000005 , 1000.8000000000000 ]
z: [ -999.62318367346938 , 1000.3768163265306 ]
If the system is larger than this, increase the sidemax parameter.
################################################################################
Building initial approximation ...
################################################################################
Adjusting initial point to fit the constraints
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Molecules of type: 1
Packing:|0 100%|
|*********************
Restraint-only function value: 1.6124008673254090E-002
Maximum violation of the restraints: 2.1233568887108977E-003
--------------------------------------------------------------------------------
Molecules of type: 2
Packing:|0 100%|
|*************************************************
Restraint-only function value: 3.7048417251638570E-005
Maximum violation of the restraints: 3.5406668572434895E-005
--------------------------------------------------------------------------------
Rescaling maximum and minimum coordinates...
Computing size of patches...
Add fixed molecules to permanent arrays...
Reseting center of mass...
--------------------------------------------------------------------------------
Setting initial trial coordinates ...
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Molecules of type: 1
Adjusting random positions to fit the constraints.
Packing:|0 100%|
|*******
Restraint-only function value: 1.6474012272086560E-006
Maximum violation of the restraints: 1.6474012272086560E-006
--------------------------------------------------------------------------------
Molecules of type: 2
Adjusting random positions to fit the constraints.
Packing:|0 100%|
|*******
Restraint-only function value: 0.0000000000000000
Maximum violation of the restraints: 0.0000000000000000
################################################################################
Objective function at initial point: 56739.880081549039
################################################################################
Packing molecules of type: 1
################################################################################
--------------------------------------------------------------------------------
Starting GENCAN loop: 0
Scaling radii by: 1.1000000000000001
Packing:|0 100%|
|*************************************************
Function value from last GENCAN loop: f = .11026E-05
Best function value before: f = .33162E+04
Improvement from best function value: 99.99 %
Improvement from last loop: 99.99 %
Maximum violation of target distance: 0.000000
Maximum violation of the constraints: .11026E-05
Current structure written to file: /home/runner/work/MDMCv0.2_pilot/MDMCv0.2_pilot/doc/how-to/use-MDMC/notebooks/packmol_files/output-universe.pdb
--------------------------------------------------------------------------------
Packing solved for molecules of type 1
Objective function value: 1.1026406870367777E-006
Maximum violation of target distance: 0.0000000000000000
Max. constraint violation: 1.1026406870367777E-006
--------------------------------------------------------------------------------
################################################################################
Packing molecules of type: 2
################################################################################
--------------------------------------------------------------------------------
Starting GENCAN loop: 0
Scaling radii by: 1.1000000000000001
Packing:|0 100%|
|******************************************************************|
|******************************************************************|
Function value from last GENCAN loop: f = .27024E+02
Best function value before: f = .33509E+05
Improvement from best function value: 99.92 %
Improvement from last loop: 99.92 %
Maximum violation of target distance: 0.615034
Maximum violation of the constraints: .14727E+01
All-type function value: .20513E+05
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Starting GENCAN loop: 1
Scaling radii by: 1.1000000000000001
Packing:|0 100%|
|******************************************************************|
|******************************************************************|
Function value from last GENCAN loop: f = .13447E+00
Best function value before: f = .27024E+02
Improvement from best function value: 99.50 %
Improvement from last loop: 99.50 %
Maximum violation of target distance: 0.000000
Maximum violation of the constraints: .45791E-01
All-type function value: .20818E+05
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Starting GENCAN loop: 2
Scaling radii by: 1.1000000000000001
Packing:|0 100%|
|********************************************************
Function value from last GENCAN loop: f = .20567E-01
Best function value before: f = .13447E+00
Improvement from best function value: 84.71 %
Improvement from last loop: 84.71 %
Maximum violation of target distance: 0.000000
Maximum violation of the constraints: .90353E-02
Current structure written to file: /home/runner/work/MDMCv0.2_pilot/MDMCv0.2_pilot/doc/how-to/use-MDMC/notebooks/packmol_files/output-universe.pdb
--------------------------------------------------------------------------------
Packing solved for molecules of type 2
Objective function value: 2.0566542813648576E-002
Maximum violation of target distance: 0.0000000000000000
Max. constraint violation: 9.0352908336764973E-003
--------------------------------------------------------------------------------
################################################################################
Packing all molecules together
################################################################################
--------------------------------------------------------------------------------
Starting GENCAN loop: 0
Scaling radii by: 1.1000000000000001
Packing:|0 100%|
|******************************************************************|
|******************************************************************|
Function value from last GENCAN loop: f = .14167E+03
Best function value before: f = .20825E+05
Improvement from best function value: 99.32 %
Improvement from last loop: 99.32 %
Maximum violation of target distance: 2.762196
Maximum violation of the constraints: .30606E+01
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Starting GENCAN loop: 1
Scaling radii by: 1.1000000000000001
Packing:|0 100%|
|******************************************************************|
|******************************************************************|
Function value from last GENCAN loop: f = .21075E+02
Best function value before: f = .14167E+03
Improvement from best function value: 85.12 %
Improvement from last loop: 85.12 %
Maximum violation of target distance: 1.251296
Maximum violation of the constraints: .10843E+01
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Starting GENCAN loop: 2
Scaling radii by: 1.1000000000000001
Packing:|0 100%|
|******************************************************************|
|******************************************************************|
Function value from last GENCAN loop: f = .27692E+01
Best function value before: f = .21075E+02
Improvement from best function value: 86.86 %
Improvement from last loop: 86.86 %
Maximum violation of target distance: 0.000000
Maximum violation of the constraints: .24400E+00
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Starting GENCAN loop: 3
Scaling radii by: 1.1000000000000001
Packing:|0 100%|
|******************************************************************|
|******************************************************************|
Function value from last GENCAN loop: f = .30068E+00
Best function value before: f = .27692E+01
Improvement from best function value: 89.14 %
Improvement from last loop: 89.14 %
Maximum violation of target distance: 0.000000
Maximum violation of the constraints: .61792E-01
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Starting GENCAN loop: 4
Scaling radii by: 1.1000000000000001
Packing:|0 100%|
|******************************************************************|
|******************************************************************|
Function value from last GENCAN loop: f = .48890E-01
Best function value before: f = .30068E+00
Improvement from best function value: 83.74 %
Improvement from last loop: 83.74 %
Maximum violation of target distance: 0.000000
Maximum violation of the constraints: .14514E-01
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Starting GENCAN loop: 5
Scaling radii by: 1.1000000000000001
Packing:|0 100%|
|***********************************
Function value from last GENCAN loop: f = .30659E-01
Best function value before: f = .48890E-01
Improvement from best function value: 37.29 %
Improvement from last loop: 37.29 %
Maximum violation of target distance: 0.000000
Maximum violation of the constraints: .94141E-02
################################################################################
Success!
Final objective function value: .30659E-01
Maximum violation of target distance: 0.000000
Maximum violation of the constraints: .94141E-02
--------------------------------------------------------------------------------
Please cite this work if Packmol was useful:
L. Martinez, R. Andrade, E. G. Birgin, J. M. Martinez,
PACKMOL: A package for building initial configurations for
molecular dynamics simulations.
Journal of Computational Chemistry, 30:2157-2164,2009.
################################################################################
Solution written to file: /home/runner/work/MDMCv0.2_pilot/MDMCv0.2_pilot/doc/how-to/use-MDMC/notebooks/packmol_files/output-universe.pdb
--------------------------------------------------------------------------------
Running time: 7.66204500 seconds.
--------------------------------------------------------------------------------
Now we can use view to look at the solved system. This may be slow - switch the viewer to ASE if it is not loading.
[5]:
from MDMC.gui import view
view(universe)
[5]:
Creating bilayers
PackmolSetup boxes, cubes, and spheres can be given a centre, for objects not centred on the origin.
Here we will demonstrate how this can be used to create an interface between benzene and water. We get our benzene molecule from file.
[6]:
from MDMC.readers.configurations import read
setup_2 = PackmolSetup()
benzene_atoms = read("data/benzene.pdb")
benzene_molecule = Molecule(atoms=benzene_atoms, name='benzene')
We define the centre of our boxes through the origin parameter;
[7]:
setup_2.add_box(structure=water_molecule, lengths=(20., 10., 20.,), origin=(0., 10., 0.), density=0.01)
setup_2.add_box(structure=benzene_molecule, lengths=(20., 10., 20.,), origin=(0., 20., 0.), density=0.01)
and as before, we fill with the PackmolFiller.
[8]:
filler_2 = PackmolFiller(setup_data=setup_2)
universe_2 = filler_2.fill_with_packmol()
################################################################################
PACKMOL - Packing optimization for the automated generation of
starting configurations for molecular dynamics simulations.
Version 20.14.2
################################################################################
Packmol must be run with: packmol < inputfile.inp
Userguide at: http://m3g.iqm.unicamp.br/packmol
Reading input file... (Control-C aborts)
Types of coordinate files specified: pdb
Seed for random number generator: 1234567
Output file: /home/runner/work/MDMCv0.2_pilot/MDMCv0.2_pilot/doc/how-to/use-MDMC/notebooks/packmol_files/output-universe.pdb
Reading coordinate file: water-0.pdb
Reading coordinate file: benzene-1.pdb
Number of independent structures: 2
The structures are:
Structure 1 :water-0.pdb( 3 atoms)
Structure 2 :benzene-1.pdb( 12 atoms)
Maximum number of GENCAN loops for all molecule packing: 400
Total number of restrictions: 2
Distance tolerance: 2.0000000000000000
Warning: Type of residue numbering not set for structure 1
Residue numbering set for structure 1 : 0
Swap chains of molecules of structure 1 : F
Warning: Type of residue numbering not set for structure 2
Residue numbering set for structure 2 : 0
Swap chains of molecules of structure 2 : F
Number of molecules of type 1 : 40
Number of molecules of type 2 : 40
Total number of atoms: 600
Total number of molecules: 80
Number of fixed molecules: 0
Number of free molecules: 80
Number of variables: 480
Total number of fixed atoms: 0
Maximum internal distance of type 1 : 1.6330000000000000
Maximum internal distance of type 2 : 4.9950000000000001
All atoms must be within these coordinates:
x: [ -1000.0000000000000 , 1000.0000000000000 ]
y: [ -987.20000000000005 , 1012.8000000000000 ]
z: [ -999.83587555555550 , 1000.1641244444445 ]
If the system is larger than this, increase the sidemax parameter.
################################################################################
Building initial approximation ...
################################################################################
Adjusting initial point to fit the constraints
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Molecules of type: 1
Packing:|0 100%|
|***********************************
Restraint-only function value: 7.7350458376762299E-003
Maximum violation of the restraints: 4.0061570777049161E-003
--------------------------------------------------------------------------------
Molecules of type: 2
Packing:|0 100%|
|********************************************************
Restraint-only function value: 4.3272511699121372E-003
Maximum violation of the restraints: 3.7571227894299315E-003
--------------------------------------------------------------------------------
Rescaling maximum and minimum coordinates...
Computing size of patches...
Add fixed molecules to permanent arrays...
Reseting center of mass...
--------------------------------------------------------------------------------
Setting initial trial coordinates ...
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Molecules of type: 1
Adjusting random positions to fit the constraints.
Packing:|0 100%|
|
Restraint-only function value: 2.1362347909628801E-003
Maximum violation of the restraints: 2.1362347909628801E-003
--------------------------------------------------------------------------------
Molecules of type: 2
Adjusting random positions to fit the constraints.
Packing:|0 100%|
|**************
Restraint-only function value: 0.0000000000000000
Maximum violation of the restraints: 0.0000000000000000
################################################################################
Objective function at initial point: 5545.9493971182119
################################################################################
Packing molecules of type: 1
################################################################################
--------------------------------------------------------------------------------
Starting GENCAN loop: 0
Scaling radii by: 1.1000000000000001
Packing:|0 100%|
|*******
Function value from last GENCAN loop: f = .52662E-02
Best function value before: f = .12277E+03
Improvement from best function value: 99.99 %
Improvement from last loop: 99.99 %
Maximum violation of target distance: 0.000000
Maximum violation of the constraints: .35041E-02
Current structure written to file: /home/runner/work/MDMCv0.2_pilot/MDMCv0.2_pilot/doc/how-to/use-MDMC/notebooks/packmol_files/output-universe.pdb
--------------------------------------------------------------------------------
Packing solved for molecules of type 1
Objective function value: 5.2662371886161236E-003
Maximum violation of target distance: 0.0000000000000000
Max. constraint violation: 3.5040911819184282E-003
--------------------------------------------------------------------------------
################################################################################
Packing molecules of type: 2
################################################################################
--------------------------------------------------------------------------------
Starting GENCAN loop: 0
Scaling radii by: 1.1000000000000001
Packing:|0 100%|
|******************************************************************|
|******************************************************************|
Function value from last GENCAN loop: f = .14687E+02
Best function value before: f = .54220E+04
Improvement from best function value: 99.73 %
Improvement from last loop: 99.73 %
Maximum violation of target distance: 0.168880
Maximum violation of the constraints: .15825E+01
All-type function value: .57528E+02
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Starting GENCAN loop: 1
Scaling radii by: 1.1000000000000001
Packing:|0 100%|
|******************************************************************|
|******************************************************************|
Function value from last GENCAN loop: f = .18740E+01
Best function value before: f = .14687E+02
Improvement from best function value: 87.24 %
Improvement from last loop: 87.24 %
Maximum violation of target distance: 0.000000
Maximum violation of the constraints: .37438E+00
All-type function value: .28198E+02
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Starting GENCAN loop: 2
Scaling radii by: 1.1000000000000001
Packing:|0 100%|
|******************************************************************|
|******************************************************************|
Function value from last GENCAN loop: f = .44152E+00
Best function value before: f = .18740E+01
Improvement from best function value: 76.44 %
Improvement from last loop: 76.44 %
Maximum violation of target distance: 0.000000
Maximum violation of the constraints: .14317E+00
All-type function value: .16619E+02
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Starting GENCAN loop: 3
Scaling radii by: 1.1000000000000001
Packing:|0 100%|
|******************************************************************|
|******************************************************************|
Function value from last GENCAN loop: f = .67993E-01
Best function value before: f = .44152E+00
Improvement from best function value: 84.60 %
Improvement from last loop: 84.60 %
Maximum violation of target distance: 0.000000
Maximum violation of the constraints: .17759E-01
All-type function value: .17511E+02
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Starting GENCAN loop: 4
Scaling radii by: 1.1000000000000001
Packing:|0 100%|
|******************************************
Function value from last GENCAN loop: f = .27235E-01
Best function value before: f = .67993E-01
Improvement from best function value: 59.94 %
Improvement from last loop: 59.94 %
Maximum violation of target distance: 0.000000
Maximum violation of the constraints: .87280E-02
Current structure written to file: /home/runner/work/MDMCv0.2_pilot/MDMCv0.2_pilot/doc/how-to/use-MDMC/notebooks/packmol_files/output-universe.pdb
--------------------------------------------------------------------------------
Packing solved for molecules of type 2
Objective function value: 2.7235229699771660E-002
Maximum violation of target distance: 0.0000000000000000
Max. constraint violation: 8.7280214800163659E-003
--------------------------------------------------------------------------------
################################################################################
Packing all molecules together
################################################################################
--------------------------------------------------------------------------------
Starting GENCAN loop: 0
Scaling radii by: 1.1000000000000001
Packing:|0 100%|
|******************************************************************|
|
Function value from last GENCAN loop: f = .31605E-01
Best function value before: f = .18006E+02
Improvement from best function value: 99.82 %
Improvement from last loop: 99.82 %
Maximum violation of target distance: 0.000000
Maximum violation of the constraints: .72331E-02
################################################################################
Success!
Final objective function value: .31605E-01
Maximum violation of target distance: 0.000000
Maximum violation of the constraints: .72331E-02
--------------------------------------------------------------------------------
Please cite this work if Packmol was useful:
L. Martinez, R. Andrade, E. G. Birgin, J. M. Martinez,
PACKMOL: A package for building initial configurations for
molecular dynamics simulations.
Journal of Computational Chemistry, 30:2157-2164,2009.
################################################################################
Solution written to file: /home/runner/work/MDMCv0.2_pilot/MDMCv0.2_pilot/doc/how-to/use-MDMC/notebooks/packmol_files/output-universe.pdb
--------------------------------------------------------------------------------
Running time: 0.246478006 seconds.
--------------------------------------------------------------------------------
[9]:
view(universe_2)
[9]: