11_Example_Input_Scripts

11 Example Input Scripts

This section contains a variety of example Ambuild input scripts, each of which involves the growth of a polymer from methane building blocks. There is one example of each of the growBlocks, joinBlocks, and zipBlocks scripts, in which the polymer is built solely from methane fragments, with no solvent, catalyst, or influence from other reagents (Sections 11.1-11.3, respectively). Additionally, some more advanced scripts are included to generate methane polymers in a statistical manner (Section 11.4) and an analysis script aimed to generate a trajectory from an Ambuild/HOOMD-blue optimisation and MD simulation that can be visualised as a DL_POLY version C HISTORY file using VMD (Section 11.5). Section 11.6 gives an example cluster generation script for Py-CMP.

11.1 Example growBlocks Input Script

#!/usr/bin/env python3

#Our imports
from ambuild import ab_cell

#Cell Dimensions
boxDim=[30,30,30]

#Create the Cell
paramsDir='./params'
mycell = ab_cell.Cell(boxDim, atomMargin=0.1, bondMargin=0.5, bondAngleMargin=5,  paramsDir=paramsDir )

#Add Fragments to the Library
mycell.libraryAddFragment( filename='./blocks/ch4.car', fragmentType='Me' )

#Specify Bonding Rules
mycell.addBondType( 'Me:a-Me:a' )

#Seed the Cell with the Building Blocks in the Library
mycell.seed( 1, fragmentType='Me' )

for i in range(100):
    mycell.growBlocks(4, cellEndGroups='Me:a', libraryEndGroups='Me:a', maxTries=500)
    mycell.optimiseGeometry(quiet=True)
    mycell.dump()
    mycell.runMD(doDihedral=True, rCut=10, mdCycles=1000000, T=55.0)

11.2 Example joinBlocks Input Script

#!/usr/bin/env python3

#Our imports
from ambuild import ab_cell

#Cell Dimensions
boxDim=[30,30,30]

#Create the Cell
paramsDir='./params'
mycell = ab_cell.Cell(boxDim, atomMargin=0.1, bondMargin=0.5, bondAngleMargin=5,  paramsDir=paramsDir )

#Add Fragments to the Library
mycell.libraryAddFragment( filename='./blocks/ch4.car', fragmentType='Me' )

#Specify Bonding Rules
mycell.addBondType( 'Me:a-Me:a' )

#Seed the Cell with the Building Blocks in the Library
mycell.seed( 100, fragmentType='Me' )
mycell.optimiseGeometry(quiet=False)

for i in range(100):
    mycell.joinBlocks(1, cellEndGroups=None, maxTries=500)
    mycell.optimiseGeometry(quiet=False)
    mycell.dump()
    mycell.runMD(doDihedral=True, rCut=10, mdCycles=1000000, T=55.0)

11.3 Example zipBlocks Input Script

#!/usr/bin/env python3

#Our imports
from ambuild import ab_cell

#Cell Dimensions
boxDim=[30,30,30]

#Create the Cell
paramsDir='./params'
mycell = ab_cell.Cell(boxDim, atomMargin=0.1, bondMargin=0.5, bondAngleMargin=5,  paramsDir=paramsDir )

#Add Fragments to the Library
mycell.libraryAddFragment( filename='./blocks/ch4.car', fragmentType='Me' )

#Specify Bonding Rules
mycell.addBondType( 'Me:a-Me:a' )

#Seed the Cell with the Building Blocks in the Library
mycell.seed( 100, fragmentType='Me' )
mycell.runMD(doDihedral=True, rCut=10, mdCycles=1000000, T=55.0)
mycell.optimiseGeometry(quiet=True)

for i in range(100):
    mycell.zipBlocks(bondMargin=2.0, bondAngleMargin=30, clashCheck=False)
    mycell.optimiseGeometry(quiet=True)
    mycell.dump()
    mycell.runMD(doDihedral=True, rCut=10, mdCycles=1000000, T=55.0)

11.4 Example Statistical Analysis Input Script

#!/usr/bin/env python3
from ambuild import ab_cell
from ambuild import ab_util

#Python imports
import csv

#Cell Dimensions
boxDim=[50,50,50]

#Create the Cell
paramsDir = 'Parameters'
mycell = ab_cell.Cell(boxDim, atomMargin=0.1, bondMargin=0.5, bondAngleMargin=5, paramsDir=paramsDir)

#Add Fragments to the Library
mycell.libraryAddFragment( filename='Blocks/ch4.car', fragmentType='CH4', markBonded=True)

#Specify Bonding Rules
mycell.addBondType( 'CH4:a-CH4:a' )

# keep a list of the configs
configs=[]

#loop to create 5 independent structures, each with 10 CH4 molecules
for i in range(5):
    #Seed
    mycell.seed(1, fragmentType='CH4', center=True)
    #Grow and zip blocks
    for i in range(9):
        mycell.growBlocks(1, cellEndGroups=['CH4:a','CH4:a'], libraryEndGroups=['CH4:a','CH4:a'])
        mycell.runMD(doDihedral=True, rCut=10, mdCycles=1000000, dt=0.0001, T=55.0, integrator='nvt')
        mycell.optimiseGeometry(doDihedral=True, rCut=10, optCycles=1000000, dt=0.0001)
        if mycell.zipBlocks(bondMargin=3, bondAngleMargin=80, selfBond=True, clashCheck=False) > 0:
            mycell.cat1Paf2(['CH4','CH4'])
            mycell.optimiseGeometry(doDihedral=True, rCut=10, optCycles=1000000, dt=0.0001)
        mycell.dump()

        # add the current endGroupConfig to the list, which gives the connectivity of each CH4 building block
        # change the text in the string (e.g. 'CH4') to the fragmentType you want to write out.
        configs.append(mycell.endGroupConfig('CH4'))
        #Printing the configs as the script goes along.
        print(configs)
    #delete CH4 blocks to start the next cluster (10 pickle steps per cluster including delete step which should be empty and shouldn't write to config file!) Removing a maximum of 10 CH4 blocks as 10 were added, but in reality only one block should be removed as these are added using the grow method.
    mycell.deleteBlocks(fragmentTypes='CH4', maxFrags=10)
    mycell.dump()

# Create the csvfile with the data we've saved
with open("endGroupConfig.csv",'w') as f:
    configWriter = csv.writer(f)
    configWriter.writerow(['FileNumber','endGroupConfig'])
    configWriter.writerows(configs)

11.5 Example Trajectory Analysis Input Script

Visualise this script using VMD, specifying the filetype as DL_POLY_C HISTORY. This script is made up of two parts: the first (First_frame_header.py) creates the 'HISTORY' file and adds the first few lines that only appear in the header of the file, and the second (VMD_analysis.py) then takes sequential Ambuild .xyz files and appends to the DL_POLY_C 'HISTORY' file by unpickling each one at a time to give the CONFIG files, then creating a HISTORY file for each frame and sticking them together to get a full HISTORY.

This needs to be run on a set of .pkl.gz files which all have the SAME number of atoms.

The two scripts are linked by the Analysis.sh file, which can be run using ./Analysis.sh. All three files need to be in the same directory as the pickle steps of interest.

Analysis.sh

#!/bin/bash

python First_frame_header.py

for i in step_*.pkl.gz
do
/opt/ambuild/ambuild/ab_util.py ${i} -da
python VMD_analysis.py ${i}
done

rm CONTROL CONFIG FIELD

To use for a different system, you will need to edit the forcefield types, masses and charges in VMD_analysis.py. The number of MD cycles, integration timestep (dt), number of atoms and unit cell length need to be inputted to the VMD_analysis.py script. The number of atoms also needs to be inputted into First_frame_header.py. If the number of atoms is too low or not all of the required forcefield types are added, the script will not include all of the atoms within the HISTORY file, which will cause it to not load successfully in VMD.

It is necessary to run this on a machine that has Ambuild installed as otherwise it won't be able to generate the DL_POLY inputs. When the script is running, Ambuild may ask for additional parameters to be added to the parameters directory, but we are not actually running the DL_POLY job, just visualising the HOOMD-blue trajectory that already contains the necessary parameters, this is not required.

First_frame_header.py

#Change the number of atoms here as appropriate.
number_of_atoms = 100

with open('HISTORY', 'w') as output:
    output.write('Ambuild CONFIG file\n         0         1       '+str(number_of_atoms)+'\n')

VMD_analysis.py

import sys
import re

#Take the input filename from the terminal in the format 'step_x.pkl.gz'
input_filename = sys.argv[1]

#Change the number of MD cycles, dt (integration timestep), number of atoms and cell_length here as appropriate. We think the standard dt=0.0001. Keep the speech marks around the cell_lengths to get them to write to the correct number of decimal places.
MD_cycles = 1000
integration_timestep = 0.0001
number_of_atoms = 100
x_cell_length = '70'
y_cell_length = '70'
z_cell_length = '200'

#Split the filename up so you just keep the pickle step number x
delimiters = '_', '.'
regex_pattern = '|'.join(map(re.escape, delimiters))
step_number = re.split(regex_pattern, input_filename)[1]
timestep = int(step_number) * int(MD_cycles)

#Open the CONFIG file as read-only and read in the data. Create an empty list called filedata to add the bits we want to.
with open('CONFIG', 'r') as CONFIG_file:
    data = CONFIG_file.readlines()

filedata = []

#Skip the first five lines as we're going to put these in manually. Then, work out which lines have the forcefield types on (as we'll need to add to these) and which have the coordinates. Then, specify the masses and charges for all of the FF types and tell it to replace the original FF line of 'c' with 'c, atom number, mass, charge' and then add the positions underneath as they are. Add all this info to filedata.
#EDIT the masses and charges to be appropriate to your system, as these are currently set to the atomic mass and atomic number, respectively!
#Add or remove atom types as appropriate for your system.
line = 0
atom_number = 0
for line_contents in data:
    line = line + 1
    if line > 5:
        if (line % 2) == 0:
            label = line_contents.split()[0]
            atom_number = atom_number + 1
            if label == 'h':
                mass = '1.008'
                charge = '1'
            elif label == 'c':
                mass = '12.011'
                charge = '6'
            elif label == 'xx':
                mass = '20'
                #Making up a random charge here but this needs correcting!
                charge = '10'
            else:
                break
            new_label = str(label) + '               ' + str(atom_number) + '   ' + str(mass) + '    ' + str(charge)
            filedata.append(new_label)
        if (line % 2) == 1:
            coordinates = '  ' + line_contents.split()[0] + '  ' + line_contents.split()[1] + '  ' + line_contents.split()[2]
            filedata.append(coordinates)

#Set up the new frame header
with open('HISTORY', 'a') as output:
    output.write('timestep      '+str(timestep)+'       '+str(number_of_atoms)+'         0         1    '+str(integration_timestep)+'\n   ' + str(x_cell_length) + '       0.000       0.000\n   0.000       ' + str(y_cell_length) + '       0.000\n   0.000       0.000       ' + str(z_cell_length) + '\n')

#Add the filedata for the new frame to the HISTORY file
with open('HISTORY', 'a') as output:
    for string_tuple in filedata:
        output.write(f'{string_tuple}\n')

11.6 Example Cluster Generation for Py-CMP

#!/usr/bin/env python3

# Our imports
from ambuild import ab_cell
from ambuild import ab_util

#Specify the parameters directory
paramsDir = 'Parameters'

#Specify the cell dimensions
boxDim = [60,60,60]

#Create the cell
mycell = ab_cell.Cell(boxDim, atomMargin = 0.1, bondMargin = 0.5, bondAngleMargin = 5, paramsDir = paramsDir)

#Add building block fragments to the library
mycell.libraryAddFragment(filename = 'Blocks/Pyrene.car', fragmentType = 'Pyrene', markBonded = True)
mycell.libraryAddFragment(filename = 'Blocks/DMF.car', fragmentType = 'DMF')

#Specify bonding rules
mycell.addBondType('Pyrene:a-Pyrene:a')

#Seed the cell with one tetrabromopyrene in the centre of the cell
mycell.seed( 1, fragmentType = 'Pyrene', center = True)
#Fill the remaining cell volume with DMF, optimise and run MD
mycell.seed( 3000, fragmentType = 'DMF')
mycell.optimiseGeometry(doDihedral = True, rCut = 10, optCycles = 1000000, dt = 0.0001)
mycell.runMD(doDihedral = True, rCut = 10, mdCycles = 100000, dt=0.0001, T = 1.2)

#Set up a loop to run over 9 steps, generating a cluster with ten pyrene units in the presence of DMF
for i in range(9):
    #Remove the DMF solvent
    mycell.deleteBlocks(fragmentTypes = 'DMF', maxFrags = 3000)
    #Grow one pyrene building block onto the existing polymer fragment, and re-seed the solvent
    mycell.growBlocks(toGrow = 1, cellEndGroups = None, libraryEndGroups = 'Pyrene:a', maxTries = 500)
    mycell.seed( 3000, fragmentType = 'DMF')
    #Establish whether any new bonds can form with relaxed bonding criteria, then run optimisation and MD
    mycell.zipBlocks(bondMargin = 5.0, bondAngleMargin = 70)
    mycell.optimiseGeometry(doDihedral = True, rCut = 10, optCycles = 1000000, dt = 0.0001)
    mycell.runMD(doDihedral = True, rCut = 10, mdCycles = 100000, dt=0.0001, T = 1.2)