auto_namd_python/simulate.py

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