保證cpu至少要有24個核 或者自己改一下代碼

#!/usr/bin/env python # -*- coding: utf-8 -*-# ============================================================================= # 對于單個復合物進行md流程 # =============================================================================import os, sys, glob, re from rdkit import Chemdef write_file(output_file, outline):buffer = open(output_file, 'w')buffer.write(outline)buffer.close()def del_conect(file_path):'''Delete contect infomation in pdb file'''lines = [l for l in open(file_path, 'r') if 'CONECT' not in l]f = open(f'{file_path}_new', 'w', encoding="utf-8")f.writelines(lines)f.closeos.remove(file_path)os.rename(f'{file_path}_new', file_path)def protein_prepare(res_name):'''this script is used for protein preparation.'''del_conect(f'{res_name}.pdb')cmdline = 'pdb4amber -i %s.pdb -o %s_pre.pdb' % (res_name, res_name)# cmdline = 'pdb2pqr %s %s.pqr --ff=amber --ffout=amber --ph-calc-method=propka --with-ph=7.0' % (protein_file, protein_name)os.system(cmdline)def get_formalcharge(infile):mol = Chem.SDMolSupplier(infile) mol = Chem.AddHs(mol) FChar = [atom.GetFormalCharge() for atom in mol.GetAtoms()] print ('the formal charge of ligand is %s'%(str(sum(FChar))))return sum(FChar)def get_resp_file(name):'''add 'resp' charges to ligand'''net_charge = get_formalcharge('%s_l.sdf' % name)# cmdline = 'cd %s &&' % namecmdline = 'antechamber -i %s_l.mol2 -fi mol2 -o %s_l.gif -fo gcrt -ch "%s_l.chk" -nc %s ' % (name, name, name, net_charge)cmdline += '-gm "%mem=2048MB" -gn "%nproc=28" &&'cmdline += 'g16 %s_l.gif > %s_l.log &&' % (name, name)cmdline += 'antechamber -i %s_l.log -fi gout -o %s_l_pre.mol2 -fo mol2 -c resp -pf y -rn "MOL" -nc %s &&' % (name, name, net_charge)cmdline += 'antechamber -i %s_l.mol2 -fi mol2 -o %s_l_ant.mol2 -fo mol2 -pf y -rn "MOL" -nc %s &&' % (name, name, net_charge)cmdline += 'antechamber -i %s_l_ant.mol2 -fi mol2 -a %s_l_pre.mol2 -fa mol2 -ao crg -o %s_l_new.mol2 -fo mol2 -pf y &&' % (name, name, name)cmdline += 'parmchk2 -i %s_l_new.mol2 -f mol2 -o %s_l_pre.frcmod' % (name, name)os.system(cmdline)def tleap_afterresp(name):#outline = 'source oldff/leaprc.ff14SB\n'#outline = 'source leaprc.protein.ff14SB\n'outline = 'source leaprc.protein.ff03.r1\n'outline += 'source leaprc.gaff\n'outline += 'loadamberparams %s_l_pre.frcmod\n' % nameoutline += 'ligand = loadmol2 %s_l_new.mol2\n' % nameoutline += 'saveoff ligand %s_l.lib\n' % nameoutline += 'saveamberparm ligand %s_l.top %s_l.crd\n' % (name, name)outline += 'quit'write_file('leap_ligand.in', outline)##outline2 = 'source oldff/leaprc.ff14SB\n'#outline2 = 'source leaprc.protein.ff14SB\n'outline2 = 'source leaprc.protein.ff03.r1\n'outline2 += 'source leaprc.gaff\n'##outline2 += 'loadamberparams frcmod.ionsjc_tip3p\n'outline2 += 'source leaprc.water.tip3p\n'#outline2 += 'receptor1 = loadpdb %s_p_pre.pdb\n' % nameoutline2 += 'receptor1 = loadpdb %s_p.pqr\n' % nameoutline2 += 'saveamberparm receptor1 %s_p.top %s_p.crd\n' % (name, name)outline2 += 'loadamberparams %s_l_pre.frcmod\n' % nameoutline2 += 'loadoff %s_l.lib\n' % nameoutline2 += 'complex = combine {receptor1 ligand}\n'outline2 += 'saveamberparm complex %s_c.top %s_c.crd\n' % (name, name)outline2 += 'addions complex Cl- 0\n'outline2 += 'addions complex Na+ 0\n'outline2 += 'charge complex\n'outline2 += 'set default PBradii mbondi2\n'outline2 += 'solvatebox complex TIP3PBOX 12.0\n'outline2 += 'saveamberparm complex %s_sol.top %s_sol.crd\n' % (name, name)outline2 += 'savepdb complex %s_sol_leap.pdb\n' % nameoutline2 += 'quit'write_file('leap_sol.in', outline2)# cmdline = 'cd %s &&' % namecmdline = 'tleap -f leap_ligand.in &&'cmdline += 'tleap -f leap_sol.in'os.system(cmdline)def get_res_num(name):'''get total residue number from leap file'''filelines = open('%s_sol_leap.pdb' % name, 'r').readlines()i = 0while filelines[i].strip() != 'TER':i = i + 1res_num = filelines[i + 1].split()[4]return res_numdef minimization(name):'''perform minimization'''res_total_num = get_res_num(name)outline1 = '''01_min_restrain.in: minimization with Cartesian restraints&cntrlimin = 1, maxcyc = 5000, ncyc = 2500,cut = 10.0, ntb = 1ntc = 2, ntf = 2,ntpr = 100,ntr = 1, restraintmask = ':1-%s&!@H=', restraint_wt = 5.0/ '''%res_total_numwrite_file('01_min_solv.in', outline1)outline2 = '''02_min_restrain.in: minimization with Cartesian restraints&cntrlimin = 1, maxcyc = 5000, ncyc = 2500,cut = 10.0, ntb = 1ntc = 2, ntf = 2,ntpr = 100,ntr = 1, restraintmask = ':1-%s@CA | :%s', restraint_wt = 4.0/ '''%(int(res_total_num)-1, res_total_num) write_file('02_min_sidechain.in', outline2)outline3 = '''03_min_restrain.in: minimization with Cartesian restraints&cntrlimin = 1, maxcyc = 5000, ncyc = 2500,cut = 10.0, ntb = 1ntc = 2, ntf = 2,ntpr = 100,/ ''' write_file('03_min_all.in', outline3) # cmdline = 'cd %s &&' % namecmdline = 'mpirun -np 28 pmemd.MPI -O -i 01_min_solv.in -o 01_min_solv.out -p %s_sol.top -c %s_sol.crd -r 01_min_solv.rst -ref %s_sol.crd -inf 01_min_solv.mdinfo -x 01_min_solv.nc &&' % (name, name, name)cmdline += 'mpirun -np 28 pmemd.MPI -O -i 02_min_sidechain.in -o 02_min_sidechain.out -p %s_sol.top -c 01_min_solv.rst -r 02_min_sidechain.rst -ref 01_min_solv.rst -inf 02_min_sidechain.mdinfo -x 02_min_sidechain.nc &&' % namecmdline += 'mpirun -np 28 pmemd.MPI -O -i 03_min_all.in -o 03_min_all.out -p %s_sol.top -c 02_min_sidechain.rst -r 03_min_all.rst -ref 02_min_sidechain.rst -inf 03_min_all.mdinfo -x 03_min_all.nc &&' % namecmdline += 'ambpdb -p %s_sol.top -c 03_min_all.rst > %s_after_min.pdb ' % (name, name)os.system(cmdline) def heat(name):res_total_num = get_res_num(name)outline = '''04_heat_100ps.in: 100ps of heat equilibration&cntrlimin = 0, irest = 0, ntx = 1,nstlim = 50000, dt = 0.002,ntc = 2, ntf = 2,cut = 8.0, ntb = 1,ntpr = 2000, ntwx = 2000,ntt = 3, gamma_ln = 2.0,tempi = 0.0, temp0 = 300.0,ntr = 1, restraintmask = ':1-%s&!@H=', restraint_wt = 3.0,nmropt = 1,iwrap = 1, ioutfm = 1,/&wt TYPE='TEMP0', istep1=0, istep2=25000,value1=0.1, value2=300.0, /&wt TYPE='END' / '''%res_total_num write_file('04_heat_100ps.in', outline) # cmdline = 'cd %s &&' % namecmdline = 'mpirun -np 28 pmemd.MPI -O -i 04_heat_100ps.in -o 04_heat_100ps.out -p %s_sol.top -c 03_min_all.rst -r 04_heat_100ps.rst -ref 03_min_all.rst -inf 04_heat_100ps.mdinfo -x 04_heat_100ps.nc ' % nameos.system(cmdline) def density(name):res_total_num = get_res_num(name)outline = '''05_density_100ps.in: 100 ps of density equilibration&cntrlimin = 0, irest = 1, ntx = 5,nstlim = 50000, dt = 0.002,ntc = 2, ntf = 2,cut = 8.0, ntb = 2, ntp = 1, taup = 1.0,ntpr = 2000, ntwx = 2000,ntt = 3, gamma_ln = 2.0,temp0 = 300.0,ntr = 1, restraintmask = ':1-%s&!@H=', restraint_wt = 2.0,iwrap = 1, ioutfm=1,/ '''%res_total_num write_file('05_density_100ps.in', outline) # cmdline = 'cd %s &&' % namecmdline = 'mpirun -np 28 pmemd.MPI -O -i 05_density_100ps.in -o 05_density_100ps.out -p %s_sol.top -c 04_heat_100ps.rst -r 05_density_100ps.rst -ref 04_heat_100ps.rst -inf 05_density_100ps.mdinfo -x 05_density_100ps.nc ' % nameos.system(cmdline)def equilibrium(name):res_total_num = get_res_num(name)outline = '''06_equil_100ps.in: 100ps of constant pressure equilibration at 300K&cntrlimin = 0, irest = 1, ntx = 5,nstlim = 50000, dt = 0.002,ntc = 2, ntf = 2,cut = 8.0, ntb = 2, ntp = 1, taup = 2.0,ntpr = 2000, ntwx = 2000,ntt = 3, gamma_ln = 2.0,temp0 = 300.0,iwrap = 1, ioutfm=1,/ '''write_file('06_equil_100ps.in', outline) cmdline = 'mpirun -np 28 pmemd.MPI -O -i 06_equil_100ps.in -o 06_equil_100ps.out -p %s_sol.top -c 05_density_100ps.rst -r 06_equil_100ps.rst -ref 05_density_100ps.rst -inf 06_equil_100ps.mdinfo -x 06_equil_100ps.nc ' % nameos.system(cmdline)def prod(name):outline = '''07_prod.in: 20ns of md production&cntrlimin = 0, irest = 1, ntx = 5,nstlim = 10000000, dt = 0.002,ntc = 2, ntf = 2,cut = 8.0, ntb = 2, ntp = 1, taup = 2.0,ntpr = 5000, ntwx = 5000, ntwr=10000000,ntt = 3, gamma_ln = 2.0, ig = -1,temp0 = 300.0,iwrap = 1, ioutfm = 1,/ '''write_file('07_prod.in', outline) cmdline = 'export CUDA_VISIBLE_DEVICES=1 &&'cmdline += 'pmemd.cuda -O -i 07_prod.in -o 07_prod.out -p %s_sol.top -c 06_equil_100ps.rst -r 07_prod.rst -ref 06_equil_100ps.rst -inf 07_prod.mdinfo -x 07_prod.nc ' % nameos.system(cmdline)def main(protein_file, mol_file):##convert.py#protein_prepare(protein_file)name = os.path.basename(mol_file).split('_')[0]#get_resp_file(name)tleap_afterresp(name)minimization(name)heat(name)density(name)equilibrium(name)#prod(name)if __name__ == '__main__':protein_file = sys.argv[1]mol_file = sys.argv[2]main(protein_file, mol_file)

總結

以上是生活随笔為你收集整理的AMBER:对单个复合物进行分子动力学模拟的python包(resp计算电荷及gpu加速版本)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。