Initial commit

This commit is contained in:
cryobry
2019-07-06 21:13:14 -04:00
commit 60ef1111c6
6 changed files with 509 additions and 0 deletions

4
.gitignore vendored Normal file
View File

@@ -0,0 +1,4 @@
auto_namd/
__pycache__
testing
hosts.sh

28
functions.py Normal file
View File

@@ -0,0 +1,28 @@
#!/usr/bin/env python3
import os
import re
import argparse
def natural_key(string_):
return [int(s) if s.isdigit() else s for s in re.split(r'(\d+)', string_)]
def abs_path(path):
abs_path = os.path.abspath(path)
return abs_path
def parser():
parser = argparse.ArgumentParser()
parser.add_argument('--steps', type=int,
help='Number of steps to run the simulation')
parser.add_argument('--jobs_path', type=str, nargs='+', action='append',
help='Jobs path containing the jobs to be run')
parser.add_argument('--namd_params', type=str,
help='Optional parameters to send to namd')
parser.add_argument('--ffs_path', type=str,
help='Location of the forcefield files')
parser.add_argument('--namd2_bin', type=str,
help='Location of the namd2 executable')
args = parser.parse_args()
return args

152
job.py Normal file
View File

@@ -0,0 +1,152 @@
import os
import sys
import glob
from functions import natural_key, abs_path
class Job:
"""
A class that represents a job/system residing in a unique directory path
Calling prepare_sim() assigns the additional instance variables:
steps = number of steps to run the simulation
ffs = a list of forcefield files for the simulation
next_step, final_step_num = the type of step and the final number of steps
conf, out file = the configuration file and the simulation output file
"""
def __init__(self, path):
self.set_empty_vars()
self.path = path
self.name = self.get_job_name()
self.prefix = os.path.join(self.path, self.name)
self.pdb, self.psf = self.get_pdb_and_psf()
self.stage, self.step, self.coor = self.get_stage()
def set_empty_vars(self):
self.pdb, self.psf = '', ''
self.steps = ''
self.ffs_path, self.ffs = '', ''
self.next_stage, self.next_step = '', ''
self.conf, self.out = '', ''
def get_job_name(self):
'''Get job name from path'''
name = os.path.basename(self.path)
return name
def get_pdb_and_psf(self):
pdb = glob.glob(f'{self.prefix}*_solv_ion.pdb')
psf = glob.glob(f'{self.prefix}*_solv_ion.psf')
if len(pdb) == 1 and len(psf) == 1:
return pdb[0], psf[0]
pdb = glob.glob(f'{self.prefix}.pdb')
psf = glob.glob(f'{self.prefix}.psf')
if len(pdb) == 1 and len(psf) == 1:
return pdb[0], psf[0]
else:
print('No PDB or PSF files found, exiting...')
sys.exit(1)
def get_stage(self):
coors_path = self.prefix + '*.coor'
coors = glob.glob(coors_path)
coors.sort(key=natural_key)
if len(coors) >= 1:
coor = coors[-1]
stage = coor.split('_')[-2].split('.')[0]
step = int(coor.split('_')[-1].split('.')[0])
else:
coor = ''
stage = ''
step = 0
return stage, int(step), coor
def prepare_sim(self, ffs_path, steps):
self.steps = steps
self.ffs_path = abs_path(ffs_path)
self.ffs = self.get_ffs()
self.next_stage, self.next_step = self.get_next_stage()
self.conf, self.out = self.conf_out_file()
def get_next_stage(self):
if self.stage == '4-sim':
next_stage = '4-sim'
next_step = self.step + self.steps
elif self.stage == '3-heat':
next_stage = '4-sim'
next_step = steps
elif self.stage == '2-min':
next_stage = '3-heat'
next_step = 0
elif self.stage == '1-min':
next_stage = '2-min'
next_step = 0
elif self.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.ffs_path, '*'))
return ffs
def conf_out_file(self):
if self.next_stage == '4-sim':
conf = f'{self.prefix}_{self.next_stage}' \
f'_{str(self.next_step)}.conf'
out = f'{self.prefix}_{self.next_stage}' \
f'_{str(self.next_step)}.out'
else:
conf = f'{self.prefix}_{self.next_stage}.conf'
out = f'{self.prefix}_{self.next_stage}.out'
return conf, out
def info(self):
print(f'Job Path: {self.path}\n'
f'Name: {self.name}\n'
f'Working PDB: {self.pdb}\n'
f'Working PSF: {self.psf}\n'
f'Previous Step: {self.stage}\n'
f'Previous Step Number: {str(self.step)}\n'
f'Next Stage: {self.next_stage}\n'
f'Next Step Number: {str(self.next_step)}\n'
f'FF Path: {self.ffs_path}\n'
f'FFs: {self.ffs}\n'
)
def get_jobs_from_path(jobs_path):
jobs_path += '/**/*.pdb'
job_dirs = glob.glob(jobs_path, recursive=True)
job_dirs_l = []
for job_dir in job_dirs:
job_dirs_l.append(os.path.dirname(job_dir))
if len(job_dirs_l) < 1:
print('No valid jobs found in jobs path!')
sys.exit(1)
return job_dirs_l
def get_next_job(jobs):
jobs.sort(key=lambda x: (x.stage, x.step))
return jobs[0]
def create_job_instances(job_paths):
job_instances = []
for job_path in job_paths:
job_instances.append(Job(job_path))
return job_instances

3
requirements.txt Normal file
View File

@@ -0,0 +1,3 @@
git+https://github.com/Eigenstate/vmd-python.git
numpy

40
run.py Executable file
View File

@@ -0,0 +1,40 @@
#!/usr/bin/env python3
import glob
import os
from job import get_next_job, get_jobs_from_path, create_job_instances
from simulate import Simulate
from functions import abs_path, parser
def main(jobs_path, ffs_path, steps, namdbin, params):
jobs_path = abs_path(jobs_path)
job_dirs_l = get_jobs_from_path(jobs_path)
jobs = create_job_instances(job_dirs_l)
job = get_next_job(jobs)
job.prepare_sim(ffs_path, steps)
job.info()
Simulate(job, namdbin, params)
if __name__ == '__main__':
args = parser()
if args.jobs_path is not None and os.path.isdir(args.jobs_path):
jobs_path = args.jobs_path
else:
jobs_path = '/home/bryan/MD/CFTR/4-jobs'
if args.namd_params is not None:
params = args.namd_params
else:
params = '+p7'
main(jobs_path=jobs_path,
ffs_path='/home/bryan/MD/CFTR/0-forcefields',
steps=1000000,
namdbin='/home/bryan/bin/namd2',
params=params)

282
simulate.py Normal file
View File

@@ -0,0 +1,282 @@
import subprocess as sp
import vmd
class Simulate:
def __init__(self, job, namdbin, params):
self.prepare(job)
self.run_namd(job, namdbin, params)
def prepare(self, job):
if job.next_stage == '1-min':
self.min1(job)
if job.next_stage == '2-min':
self.min2(job)
if job.next_stage == '3-heat':
self.heat(job)
if job.next_stage == '4-sim':
self.sim(job)
def run_namd(self, job, namdbin, params):
cmd = f'{namdbin} {params} {job.conf} > {job.out}'
print(f'Running: {cmd}')
p1 = sp.run(cmd, shell=True)
def min1(self, job):
steps = 15000
with open(job.conf, 'w') as conf:
conf.write('# Input\n')
conf.write(f'structure {job.psf}\n')
conf.write(f'coordinates {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 {job.name}_{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, job):
#self.solv_ion(job)
cbvx, cbvy, cbvz, corx, cory, corz = self.calc_pcell(job)
steps = 15000
with open(job.conf, 'w') as conf:
conf.write('# Input\n')
conf.write(f'structure {job.psf}\n')
conf.write(f'coordinates {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('# 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 {job.name}_{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, job):
temp_reinit_steps = 100
steps = 10000
final_steps = 30 * temp_reinit_steps + steps
cbvx, cbvy, cbvz, corx, cory, corz = self.calc_pcell(job)
with open(job.conf, 'w') as conf:
conf.write('# Input\n')
conf.write(f'structure {job.psf}\n')
conf.write(f'coordinates {job.pdb}\n')
conf.write(f'bincoordinates {job.coor}\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('# 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 {job.name}_{job.next_stage}_{final_steps}\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, job):
cbvx, cbvy, cbvz, corx, cory, corz = self.calc_pcell(job)
with open(job.conf, 'w') as conf:
conf.write('# Input\n')
conf.write(f'structure {job.psf}\n')
conf.write(f'coordinates {job.pdb}\n')
conf.write(f'bincoordinates {job.coor}\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('# 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 {job.name}_{job.next_stage}_{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 {job.steps}')
def del_all_mols(self):
for mol in vmd.molecule.listall():
vmd.molecule.delete(mol)
def solv_ion(self, job):
pdb_f = f'{job.prefix}_{job.stage}.pdb'
psf_f = f'{job.prefix}.psf'
solv_f = f'{job.prefix}_{job.stage}_solv'
solv_ion_f = f'{job.prefix}_{job.stage}_solv_ion'
molid = vmd.molecule.load('psf', psf_f, 'namdbin', job.coor)
vmd.molecule.write(molid, 'pdb', pdb_f)
vmd.evaltcl('package require solvate')
vmd.evaltcl(f'solvate {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, job):
if job.next_stage == '2-min' or job.next_stage == '3-min':
molid = vmd.molecule.load('psf', job.psf, 'pdb', job.pdb)
else:
molid = vmd.molecule.load('psf', job.psf, 'namdbin', 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