123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356 |
- import glob
- import subprocess as sp
- import vmd
- from auto_namd.functions import *
- class Simulation:
- """
- A class representing a simulation of a job
- """
- def __init__(self, job, ffs_path, steps, namdbin, params):
- self.job = job
- self.steps = int(steps)
- self.namdbin = namdbin
- self.params = params
- # NAMD requires FF paths relative to the NAMD configuration file
- # and not the CWD. So we make the paths to the ffs absolute
- self.job.ffs_path = abs_path(san_path(ffs_path))
- self.job.ffs = self.get_ffs()
- self.job.next_stage, self.job.next_step = self.get_next_stage()
- self.job.conf, self.job.out = self.set_conf_out()
- print(self.info())
- self.namd_conf()
- self.namd()
- def get_next_stage(self):
- if self.job.stage == '4-sim':
- next_stage = '4-sim'
- next_step = self.job.step + self.steps
- elif self.job.stage == '3-heat':
- next_stage = '4-sim'
- next_step = self.steps
- elif self.job.stage == '2-min':
- next_stage = '3-heat'
- next_step = 0
- elif self.job.stage == '1-min':
- next_stage = '2-min'
- next_step = 0
- elif self.job.stage == '':
- next_stage = '1-min'
- next_step = 0
- return next_stage, int(next_step)
- def get_ffs(self):
- ffs = glob.glob(os.path.join(self.job.ffs_path, '*'))
- return ffs
- def set_conf_out(self):
- """Set the NAMD configration and NAMD output filenames"""
- if self.job.next_stage == '4-sim':
- conf = f'{self.job.prefix}_{self.job.next_stage}' \
- f'_{str(self.next_step)}.conf'
- out = f'{self.job.prefix}_{self.job.next_stage}' \
- f'_{str(self.job.next_step)}.out'
- else:
- conf = f'{self.job.prefix}_{self.job.next_stage}.conf'
- out = f'{self.job.prefix}_{self.job.next_stage}.out'
- return conf, out
- def namd_conf(self):
- if self.job.next_stage == '1-min':
- self.min1()
- if self.job.next_stage == '2-min':
- self.min2()
- if self.job.next_stage == '3-heat':
- self.heat()
- if self.job.next_stage == '4-sim':
- self.sim()
- def min1(self):
- steps = 15000
- with open(self.job.conf, 'w') as conf:
- conf.write('# Input\n')
- conf.write(f'structure {self.job.psf}\n')
- conf.write(f'coordinates {self.job.pdb}\n')
- conf.write('paraTypeCharmm on\n')
- for ff in job.ffs:
- conf.write(f'parameters {ff}\n')
- conf.write('\n')
- conf.write('# Temperature\n')
- conf.write('temperature 0\n')
- conf.write('\n')
- conf.write('# Force-Field Parameters\n')
- conf.write('exclude scaled1-4\n')
- conf.write('1-4scaling 1.0\n')
- conf.write('cutoff 12.\n')
- conf.write('switching on\n')
- conf.write('switchdist 10.\n')
- conf.write('pairlistdist 14\n')
- conf.write('\n')
- conf.write('# Integrator Parameters\n')
- conf.write('timestep 1.0 ;# 1fs/step\n')
- conf.write('nonbondedFreq 1\n')
- conf.write('fullElectFrequency 2\n')
- conf.write('stepspercycle 10\n')
- conf.write('\n')
- conf.write('# Output\n')
- conf.write(f'outputName {self.job.name}_{self.job.next_stage}_{steps}\n')
- conf.write('outputEnergies 100\n')
- conf.write('outputPressure 100\n')
- conf.write('\n')
- conf.write('# Run\n')
- conf.write(f'minimize {steps}')
- def min2(self):
- solv_ion(self.job)
- cbvx, cbvy, cbvz, corx, cory, corz = self.calc_pcell(self.job)
- steps = 15000
- with open(self.job.conf, 'w') as conf:
- conf.write('# Input\n')
- conf.write(f'structure {self.job.psf}\n')
- conf.write(f'coordinates {self.job.pdb}\n')
- conf.write('paraTypeCharmm on\n')
- for ff in self.job.ffs:
- conf.write(f'parameters {ff}\n')
- conf.write('\n')
- conf.write('# Temperature\n')
- conf.write('temperature 0\n')
- conf.write('\n')
- conf.write('# Periodic Boundary Conditions\n')
- conf.write('wrapWater on\n')
- conf.write('wrapAll on\n')
- conf.write(f'cellOrigin\t{corx}\t{cory}\t{corz}\n')
- conf.write(f'cellBasisVector1\t{cbvx}\t0.0\t0.0\n')
- conf.write(f'cellBasisVector2\t0.0\t{cbvy}\t0.0\n')
- conf.write(f'cellBasisVector3\t0.0\t0.0\t{cbvz}\n')
- conf.write('\n')
- conf.write('# Force-Field Parameters\n')
- conf.write('exclude scaled1-4\n')
- conf.write('1-4scaling 1.0\n')
- conf.write('cutoff 12.\n')
- conf.write('switching on\n')
- conf.write('switchdist 10.\n')
- conf.write('pairlistdist 14\n')
- conf.write('\n')
- conf.write('# Integrator Parameters\n')
- conf.write('timestep 1.0 ;# 1fs/step\n')
- conf.write('nonbondedFreq 1\n')
- conf.write('fullElectFrequency 2\n')
- conf.write('stepspercycle 10\n')
- conf.write('\n')
- conf.write('# Output\n')
- conf.write(f'outputName {self.job.name}_{self.job.next_stage}_{steps}\n')
- conf.write('outputEnergies 100\n')
- conf.write('outputPressure 100\n')
- conf.write('\n')
- conf.write('# Run\n')
- conf.write(f'minimize {steps}')
- def heat(self):
- temp_reinit_steps = 100
- steps = 10000
- final_steps = 30 * temp_reinit_steps + steps
- cbvx, cbvy, cbvz, corx, cory, corz = self.calc_pcell()
- with open(self.job.conf, 'w') as conf:
- conf.write('# Input\n')
- conf.write(f'structure {self.job.psf}\n')
- conf.write(f'coordinates {self.job.pdb}\n')
- conf.write(f'bincoordinates {self.job.coor}\n')
- conf.write('paraTypeCharmm on\n')
- for ff in self.job.ffs:
- conf.write(f'parameters {ff}\n')
- conf.write('\n')
- conf.write('# Temperature\n')
- conf.write('temperature 0\n')
- conf.write('\n')
- conf.write('# Periodic Boundary Conditions\n')
- conf.write('wrapWater on\n')
- conf.write('wrapAll on\n')
- conf.write(f'cellOrigin\t{corx}\t{cory}\t{corz}\n')
- conf.write(f'cellBasisVector1\t{cbvx}\t0.0\t0.0\n')
- conf.write(f'cellBasisVector2\t0.0\t{cbvy}\t0.0\n')
- conf.write(f'cellBasisVector3\t0.0\t0.0\t{cbvz}\n')
- conf.write('# Force-Field Parameters\n')
- conf.write('exclude scaled1-4\n')
- conf.write('1-4scaling 1.0\n')
- conf.write('cutoff 12.\n')
- conf.write('switching on\n')
- conf.write('switchdist 10.\n')
- conf.write('pairlistdist 14\n')
- conf.write('\n')
- conf.write('# Full Electrostatics\n')
- conf.write('PME on\n')
- conf.write('PMEGridSpacing 1.0\n')
- conf.write('\n')
- conf.write('# Integrator Parameters\n')
- conf.write('timestep 1.0 ;# 1fs/step\n')
- conf.write('nonbondedFreq 1\n')
- conf.write('fullElectFrequency 2\n')
- conf.write('stepspercycle 10\n')
- conf.write('\n')
- conf.write('# Output\n')
- conf.write(f'outputName {self.job.name}_{self.job.next_stage}_{self.job.next_step}\n')
- conf.write('outputEnergies 100\n')
- conf.write('outputPressure 100\n')
- conf.write('dcdfreq 1000\n')
- conf.write('\n')
- conf.write('# Constant Temperature Control\n')
- conf.write('langevin on ;# do langevin dynamics\n')
- conf.write('langevinDamping 0.5 ;# damping coefficient (gamma) of 0.5/ps\n')
- conf.write('langevinTemp 310\n')
- conf.write('langevinHydrogen yes ;# couple langevin bath to hydrogens\n')
- conf.write('\n')
- conf.write('# Constant Pressure Control\n')
- conf.write('useGroupPressure no ;# needed for 2fs steps\n')
- conf.write('useFlexibleCell yes ;# no for water box, yes for membrane\n')
- conf.write('useConstantRatio yes ;# no for water box, yes for membrane\n')
- conf.write('langevinPiston on\n')
- conf.write('langevinPistonTarget 1.01325 ;# in bar -> 1 atm\n')
- conf.write('langevinPistonPeriod 100.\n')
- conf.write('langevinPistonDecay 50.\n')
- conf.write('langevinPistonTemp 310\n')
- conf.write('\n')
- conf.write('# Run equilibration\n')
- conf.write(f'set freq {temp_reinit_steps}\n')
- conf.write('for {set i 10} {$i <= 310} {incr i 10} {\n')
- conf.write('reinitvels $i\n')
- conf.write('langevinTemp $i\n')
- conf.write('run $freq\n')
- conf.write('}\n')
- conf.write('# Run stabilization\n')
- conf.write(f'run {steps}')
- def sim(self):
- cbvx, cbvy, cbvz, corx, cory, corz = self.calc_pcell()
- with open(self.job.conf, 'w') as conf:
- conf.write('# Input\n')
- conf.write(f'structure {self.job.psf}\n')
- conf.write(f'coordinates {self.job.pdb}\n')
- conf.write(f'bincoordinates {self.job.coor}\n')
- conf.write('paraTypeCharmm on\n')
- for ff in self.job.ffs:
- conf.write(f'parameters {ff}\n')
- conf.write('\n')
- conf.write('# Temperature\n')
- conf.write('temperature 0\n')
- conf.write('\n')
- conf.write('# Periodic Boundary Conditions\n')
- conf.write('wrapWater on\n')
- conf.write('wrapAll on\n')
- conf.write(f'cellOrigin\t{corx}\t{cory}\t{corz}\n')
- conf.write(f'cellBasisVector1\t{cbvx}\t0.0\t0.0\n')
- conf.write(f'cellBasisVector2\t0.0\t{cbvy}\t0.0\n')
- conf.write(f'cellBasisVector3\t0.0\t0.0\t{cbvz}\n')
- conf.write('# Force-Field Parameters\n')
- conf.write('exclude scaled1-4\n')
- conf.write('1-4scaling 1.0\n')
- conf.write('cutoff 12.\n')
- conf.write('switching on\n')
- conf.write('switchdist 10.\n')
- conf.write('pairlistdist 14\n')
- conf.write('\n')
- conf.write('# Full Electrostatics\n')
- conf.write('PME on\n')
- conf.write('PMEGridSpacing 1.0\n')
- conf.write('\n')
- conf.write('# Integrator Parameters\n')
- conf.write('timestep 2.0 ;# 2fs/step\n')
- conf.write('rigidBonds all ;# needed for 2fs steps\n')
- conf.write('nonbondedFreq 1\n')
- conf.write('fullElectFrequency 2\n')
- conf.write('stepspercycle 10\n')
- conf.write('\n')
- conf.write('# Output\n')
- conf.write(f'outputName {self.job.name}_{self.job.next_stage}_{self.job.next_step}\n')
- conf.write('outputEnergies 10000\n')
- conf.write('outputPressure 10000\n')
- conf.write('dcdfreq 10000\n')
- conf.write('\n')
- conf.write('# Constant Temperature Control\n')
- conf.write('langevin on ;# do langevin dynamics\n')
- conf.write('langevinDamping 0.5 ;# damping coefficient (gamma) of 0.5/ps\n')
- conf.write('langevinTemp 310\n')
- conf.write('langevinHydrogen no ;# couple langevin bath to hydrogens\n')
- conf.write('\n')
- conf.write('# Constant Pressure Control\n')
- conf.write('useGroupPressure yes ;# needed for 2fs steps\n')
- conf.write('useFlexibleCell yes ;# no for water box, yes for membrane\n')
- conf.write('useConstantRatio yes ;# no for water box, yes for membrane\n')
- conf.write('langevinPiston on\n')
- conf.write('langevinPistonTarget 1.01325 ;# in bar -> 1 atm\n')
- conf.write('langevinPistonPeriod 100.\n')
- conf.write('langevinPistonDecay 50.\n')
- conf.write('langevinPistonTemp 310\n')
- conf.write('\n')
- conf.write('# Run\n')
- conf.write(f'run {self.job.steps}')
- def solv_ion(self):
- pdb_f = f'{self.job.prefix}_{self.job.stage}.pdb'
- psf_f = f'{self.job.prefix}.psf'
- solv_f = f'{self.job.prefix}_{self.job.stage}_solv'
- solv_ion_f = f'{self.job.prefix}_{self.job.stage}_solv_ion'
- molid = vmd.molecule.load('psf', psf_f, 'namdbin', self.job.coor)
- vmd.molecule.write(molid, 'pdb', pdb_f)
- vmd.evaltcl('package require solvate')
- vmd.evaltcl(f'solvate {self.job.psf} {pdb_f} -o {solv_f} '
- f'-s WT -x 13 -y 13 -z 13 +x 13 +y 13 +z 13 -b 2.4')
- vmd.evaltcl('package require autoionize')
- vmd.evaltcl(f'autoionize -psf {solv_f}.psf -pdb {solv_f}.pdb -o {solv_ion_f} '
- f'-sc 0.15')
- self.del_all_mols()
- def calc_pcell(self):
- if self.job.next_stage == '2-min' or self.job.next_stage == '3-min':
- molid = vmd.molecule.load('psf', self.job.psf, 'pdb', self.job.pdb)
- else:
- molid = vmd.molecule.load('psf', self.job.psf, 'namdbin', self.job.coor)
- all = vmd.atomsel("all", molid=molid)
- minmax = all.minmax()
- center = all.center()
- cbvx = minmax[1][0] - minmax[0][0]
- cbvy = minmax[1][1] - minmax[0][1]
- cbvz = minmax[1][2] - minmax[0][2]
- corx = center[0]
- cory = center[1]
- corz = center[2]
- self.del_all_mols()
- return cbvx, cbvy, cbvz, corx, cory, corz
- def del_all_mols(self):
- for mol in vmd.molecule.listall():
- vmd.molecule.delete(mol)
- def namd(self):
- cmd = f'{self.namdbin} {self.params} {self.job.conf} > {self.job.out}'
- print(f'Running: {cmd}')
- p1 = sp.run(cmd, shell=True)
- def info(self):
- print(f'Job Path: {self.job.path}\n'
- f'Name: {self.job.name}\n'
- f'Working PDB: {self.job.pdb}\n'
- f'Working PSF: {self.job.psf}\n'
- f'Working COOR: {self.job.coor}\n'
- f'Previous Step: {self.job.stage}\n'
- f'Previous Step Number: {str(self.job.step)}\n'
- f'Next Stage: {self.job.next_stage}\n'
- f'Next Step Number: {str(self.job.next_step)}\n'
- f'FF Path: {self.job.ffs_path}\n'
- f'FFs: {self.job.ffs}\n'
- )
|