diff --git a/bin/cmmde.py b/bin/cmmde.py index 6e1b677..91b7adc 100755 --- a/bin/cmmde.py +++ b/bin/cmmde.py @@ -58,6 +58,7 @@ parser.add_argument('-ct', '--c_terlarut', type=str, help='Muatan bersih molekul parser.add_argument('-pt', '--persen_terlarut', type=str, default=50, help='Persen massa terlarut. Jika lebih dari satu terlarut, pisahkan dengan koma.') parser.add_argument('-pp', '--persen_pelarut', type=str, default=50, help='Persen massa pelarut. Default:50') parser.add_argument('-l', '--lapang', type=float, default=5, help='Panjang tambahan (Angstrom) pada rusuk kubus untuk menghindari bad contact. Default: 10.0.') +parser.add_argument('-rc','--rcol',type=float,default=0, help='Additional Coulombic radius for PME calculation') parser.add_argument('-gen','--generate_dftbinp',type=str,default='false',help='Apakah ingin mengkonversi ke dalam format koordinat xyz?') parser.add_argument('-mp', '--pelarut', type=str, help='Nama file molekul pelarut dalam format .xyz (ditulis tanpa ekstensi).') parser.add_argument('-nequil','--nequil',type=int, default=50000, help='Jumlah step pada saat ekuilibrasi.') @@ -180,7 +181,7 @@ parser.add_argument('-cpress','--cellpress',type=float,default=0.0,help='Tekanan parser.add_argument('-press_conv','--press_conv_thr',type=float,default=0.5,help='Kriteria konvergensi untuk tekanan sel. Default=0.5 kbar.') parser.add_argument('-nband','--nband',type=int,default=8,help='Jumlah tingkat energi yang ingin diplot dalam DOS.') parser.add_argument('-occ','--occ',type=str,default='tetrahedra',help='Metode smearing yang digunakan. Default: tetrahedra. Pilihan: smearing, tetrahedra_lin, tetrahedra_opt, fixed, from_input.') -parser.add_argument("-slurm","--slurm",type=str,default="true",help="Apakah akan menggunakan slurm atau tidak. Default: true") +parser.add_argument("-q","--queue",type=str,default="true",help="Apakah akan menggunakan sistem antrian atau tidak. Default: true") opt=parser.parse_args(sys.argv[1:]) @@ -188,7 +189,7 @@ opt=parser.parse_args(sys.argv[1:]) if ('.xyz' or '.pdb' or 'POSCAR' or '.vasp' or '.poscar') in opt.input: geom = opt.input elif '.mol2' in opt.input: - if opt.slurm == "true": + if opt.queue == "true": with open('run_babel.sh','w') as fout: print("""#!/bin/bash #SBATCH --nodes=1 @@ -210,7 +211,7 @@ obabel {} -O geom.xyz""".format(opt.nproc,opt.input),file=fout) geom = 'geom.xyz' elif ('.mol2' not in opt.input and '.xyz' not in opt.input and '.pdb' not in opt.input and '.vasp' not in opt.input and '.poscar' not in opt.input and 'POSCAR' not in opt.input and '.gen' not in opt.input): os.system("echo '{}' > geom.smi".format(opt.input)) - if opt.slurm == "true": + if opt.queue == "true": with open('run_babel.sh','w') as fout: print("""#!/bin/bash #SBATCH --nodes=1 @@ -255,7 +256,7 @@ if opt.software == 'dock': addcharge(opt.input,opt.chargetype) if opt.job == 'sphgen': sphgen(opt.input) - if opt.slurm == 'true': + if opt.queue == 'true': with open('run_sphgen.sh','w') as fout: print("""#!/bin/bash #SBATCH --nodes=1 @@ -280,7 +281,7 @@ $DOCK_DIR/sphere_selector protein.sph {} {}""".format(opt.nproc,opt.ligand,opt.d showsphere() if opt.job == 'gridgen': gridgen(opt.input) - if opt.slurm == "true": + if opt.queue == "true": with open('run_grid.sh','w') as fout: print("""#!/bin/bash #SBATCH --nodes=1 @@ -302,7 +303,7 @@ $DOCK_DIR/sphere_selector protein.sph {} {}""".format(opt.nproc,opt.ligand,opt.d if opt.job == 'rigiddock': rigiddock(opt.ligand,opt.calcrmsd) - if opt.slurm == "true": + if opt.queue == "true": with open('run_rigiddock.sh','w') as fout: print("""#!/bin/bash #SBATCH --nodes=1 @@ -323,7 +324,7 @@ $DOCK_DIR/sphere_selector protein.sph {} {}""".format(opt.nproc,opt.ligand,opt.d os.system("./run_rigiddock.sh") if opt.job == 'flexdock': flexdock(opt.ligand,opt.calcrmsd) - if opt.slurm == "true": + if opt.queue == "true": with open('run_flexdock.sh','w') as fout: print("""#!/bin/bash #SBATCH --nodes=1 @@ -367,7 +368,7 @@ if opt.software == 'dftb': geom = 'in.gen' dftb(geom,opt.job,opt.activeatoms,opt.method,opt.parapath,opt.dispersion,opt.kpts,opt.hcorr) - if opt.slurm == "true": + if opt.queue == "true": with open('run.sh','w') as fout: print("""#!/bin/bash #SBATCH --nodes=1 @@ -395,7 +396,7 @@ if opt.software == 'qe': # RUNNING SCRIPT if opt.software == 'orca': - if opt.slurm=="true": + if opt.queue=="true": with open('run.sh','w') as fout: print("""#!/bin/bash #SBATCH --nodes=1 @@ -416,7 +417,7 @@ if opt.software == 'orca': os.system("chmod +x run.sh") os.system('./run.sh') if opt.software == 'nwchem': - if opt.slurm=="true": + if opt.queue=="true": with open('run.sh','w') as fout: print("""#!/bin/bash #SBATCH --nodes=1 @@ -449,7 +450,7 @@ mpirun -np {} $QE_COMMAND < cmmd.in > cmmd.out""".format(opt.nproc,opt.nproc),fi os.system('sbatch run.sh') if opt.software == 'dcdftb': - if opt.slurm == 'true': + if opt.queue == 'true': os.system("cp cmmd.in dftb.inp") with open('run.sh','w') as fout: print("""#!/bin/bash @@ -477,23 +478,23 @@ fi""",file=fout) os.system('./run.sh') if opt.software == 'gromacs': - if opt.slurm == 'true': + if opt.queue == 'true': with open('run.sh','w') as fout: print("""#!/bin/bash #SBATCH --nodes=1 #SBATCH --ntasks=1 #SBATCH --cpus-per-task=1 #SBATCH --time=168:0:0 -export OMP_NUM_THREADS={} +export OMP_NUM_THREADS=1 cd $PWD -$GROMACS_COMMAND -mt {} -mp {} -ct {} -pt {} -pp {} -l {} -cat {} -gen {} -Nump {} -prod {} -nprod {} -nequil {} -dt {} -ctype {} -nnpt {} -comp {}""".format(opt.nproc, opt.terlarut,opt.pelarut,opt.c_terlarut,opt.persen_terlarut,opt.persen_pelarut,opt.lapang,opt.cation,opt.generate_dftbinp,opt.NumPelarut,opt.production,opt.nprod,opt.nequil,opt.deltat,opt.charge_type,opt.nnpt,opt.compress),file=fout) +$GROMACS_COMMAND -mt {} -mp {} -ct {} -pt {} -pp {} -l {} -cat {} -gen {} -Nump {} -prod {} -nprod {} -nequil {} -dt {} -ctype {} -nnpt {} -comp {} -np {}""".format(opt.terlarut,opt.pelarut,opt.c_terlarut,opt.persen_terlarut,opt.persen_pelarut,opt.lapang,opt.cation,opt.generate_dftbinp,opt.NumPelarut,opt.production,opt.nprod,opt.nequil,opt.deltat,opt.charge_type,opt.nnpt,opt.compress,opt.nproc,opt.rcol),file=fout) os.system('sbatch run.sh') else: with open('run.sh','w') as fout: print("""#!/bin/bash -export OMP_NUM_THREADS={} +export OMP_NUM_THREADS=1 cd $PWD -$GROMACS_COMMAND -mt {} -mp {} -ct {} -pt {} -pp {} -l {} -cat {} -gen {} -Nump {} -prod {} -nprod {} -nequil {} -dt {} -ctype {} -nnpt {} -comp {}""".format(opt.nproc, opt.terlarut,opt.pelarut,opt.c_terlarut,opt.persen_terlarut,opt.persen_pelarut,opt.lapang,opt.cation,opt.generate_dftbinp,opt.NumPelarut,opt.production,opt.nprod,opt.nequil,opt.deltat,opt.charge_type,opt.nnpt,opt.compress),file=fout) +$GROMACS_COMMAND -mt {} -mp {} -ct {} -pt {} -pp {} -l {} -cat {} -gen {} -Nump {} -prod {} -nprod {} -nequil {} -dt {} -ctype {} -nnpt {} -comp {} -np {}""".format(opt.terlarut,opt.pelarut,opt.c_terlarut,opt.persen_terlarut,opt.persen_pelarut,opt.lapang,opt.cation,opt.generate_dftbinp,opt.NumPelarut,opt.production,opt.nprod,opt.nequil,opt.deltat,opt.charge_type,opt.nnpt,opt.compress,opt.nproc,opt.rcol),file=fout) os.system('chmod +x run.sh') os.system('./run.sh') @@ -505,4 +506,4 @@ if opt.job == 'ligprep' and opt.software == 'gmx': # Program XTB Standalone if opt.software == 'xtb': - xtb(opt.job,geom,opt.nproc,opt.produk,opt.temp,opt.nrun,opt.npoint,opt.anopt,opt.kpush,opt.kpull,opt.ppull,opt.alp,opt.distance,opt.angle,opt.dihedral,opt.scanmode,opt.iter,opt.scan,opt.solvent,opt.charge,opt.mult,opt.method,opt.fixedatoms,opt.fixedelements,opt.slurm) \ No newline at end of file + xtb(opt.job,geom,opt.nproc,opt.produk,opt.temp,opt.nrun,opt.npoint,opt.anopt,opt.kpush,opt.kpull,opt.ppull,opt.alp,opt.distance,opt.angle,opt.dihedral,opt.scanmode,opt.iter,opt.scan,opt.solvent,opt.charge,opt.mult,opt.method,opt.fixedatoms,opt.fixedelements,opt.queue) diff --git a/lib/cmmde_solution.py b/lib/cmmde_solution.py index 3858469..f4a5cc3 100755 --- a/lib/cmmde_solution.py +++ b/lib/cmmde_solution.py @@ -30,6 +30,8 @@ parser.add_argument('--packmol', type=str, default='/Users/adit/opt/packmol/pack parser.add_argument('-l','--lapang',type=float, default=5, help='Ruang lapang untuk atom-atom bergerak di dalam kotak larutan (angstrom).') parser.add_argument('-gen','--generate_dftbinp',type=str,default='false',help='Apakah ingin mengkonversi ke dalam format koordinat xyz?') parser.add_argument('-prod','--production',type=str,default='None',help='Jenis production run yang ingin dilakukan. Pilihan: NPT (recommended), NVE, dan NVT.') +parser.add_argument('-np','--nproc',type=int,default=1, help='Jumlah CPU yang digunakan') +parser.add_argument('-rc','--rcol',type=float,default=0, help='Additional Coulombic radius for PME calculation') opt=parser.parse_args(sys.argv[1:]) @@ -52,7 +54,7 @@ pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions ;define = -DFLEXIBLE """ nve_mdp = """ -define = -DFLEXIBLE +;define = -DFLEXIBLE integrator = md nsteps = {} ; 1 * 500000 = 50 ps dt = 0.001 ; 0.5 fs @@ -89,7 +91,7 @@ gen_vel = yes ; Velocity generation is off # Input file simulasi menggunakan ensembel kanonik (NVT) dalam GROMACS nvt_mdp=""" -define = -DFLEXIBLE +;define = -DFLEXIBLE integrator = md ; leap-frog integrator nsteps = {} dt = 0.001 ; 1 fs @@ -342,7 +344,7 @@ def chargeparm(mol2file,dftbdat): print(" {} {} {} {} {} {} {} {} {}".format(serial, symbol, x, y, z, type, resid, resname, charge), file=f) for line in footer: print(line, file=f) -def make_solution(terlarut,c_terlarut,pelarut,c_pelarut,N_pelarut,persen_terlarut,persen_pelarut,cation, temperature,pressure,packmol,charge_type,lapang,generate_dftbinp,prod): +def make_solution(terlarut,c_terlarut,pelarut,c_pelarut,N_pelarut,persen_terlarut,persen_pelarut,cation, temperature,pressure,packmol,charge_type,lapang,generate_dftbinp,prod,nproc): if temperature: print (f'Suhu: {temperature:12.4f} K') @@ -560,31 +562,32 @@ quit""",file=fout) nnpt = opt.nnpt nprod = opt.nprod compress = opt.compress - + rcoulomb = half_box+opt.rcol + rvdw = half_box with open('minim.mdp', 'w') as f: - print(minim_mdp.format(half_box, half_box), file=f) + print(minim_mdp.format(rcoulomb, rvdw), file=f) with open('nve.mdp', 'w') as f: - print(nve_mdp.format(nequil,half_box,half_box), file=f) + print(nve_mdp.format(nequil,rcoulomb,rvdw), file=f) with open('nvt.mdp', 'w') as f: - print(nvt_mdp.format(nequil,half_box, half_box,temperature), file=f) + print(nvt_mdp.format(nequil,rcoulomb, rvdw,temperature), file=f) with open('npt.mdp', 'w') as f: - print(npt_mdp.format(nnpt,half_box, half_box,temperature,pressure,compress), file=f) + print(npt_mdp.format(nnpt,rcoulomb, rvdw,temperature,pressure,compress), file=f) os.system("cp system.amb2gmx/system_GMX.* .") gmxscript.grompp(f="minim.mdp", c="system_box.pdb", o="mini.tpr", p="system_GMX.top", maxwarn=2) - gmxscript.mdrun(deffnm="mini",ntmpi=1) + gmxscript.mdrun(deffnm="mini",ntmpi='{}'.format(opt.nproc)) gmxscript.select(f="mini.trr", s="mini.tpr", on="mini", select="System") gmxscript.grompp(f="nve.mdp", c="mini.gro", p="system_GMX.top", o="nve.tpr",maxwarn=2) - gmxscript.mdrun(deffnm="nve", v=True, ntmpi=1) + gmxscript.mdrun(deffnm="nve", v=True, ntmpi='{}'.format(opt.nproc)) gmxscript.grompp(f="nvt.mdp", c="nve.gro", p="system_GMX.top", o="nvt.tpr",maxwarn=2) - gmxscript.mdrun(deffnm="nvt", v=True, ntmpi=1) + gmxscript.mdrun(deffnm="nvt", v=True, ntmpi='{}'.format(opt.nproc)) gmxscript.grompp(f="npt.mdp", c="nvt.gro", p="system_GMX.top", t="nvt.cpt",o="npt.tpr",maxwarn=2) - gmxscript.mdrun(deffnm="npt", v=True, ntmpi=1) + gmxscript.mdrun(deffnm="npt", v=True, ntmpi='{}'.format(opt.nproc)) if prod == 'npt': with open('nptp.mdp', 'w') as f: - print(nptp_mdp.format(nprod,opt.dt/2000,half_box,half_box,temperature,compress,pressure), file=f) + print(nptp_mdp.format(nprod,opt.dt/2000,rcoulomb,rvdw,temperature,compress,pressure), file=f) gmxscript.grompp(f="nptp.mdp", c="npt.gro", p="system_GMX.top", o="nptp.tpr",maxwarn=2) - gmxscript.mdrun(deffnm="nptp", v=True, ntmpi=1) + gmxscript.mdrun(deffnm="nptp", v=True, ntmpi='{}'.format(opt.nproc)) gmxscript.trjconv(f="nptp.trr", s="nptp.tpr", n="mini.ndx", o="finalsystem.pdb") print("Melakukan perhitungan MSD untuk pelarut") os.system("""gmx select -f nptp.gro -s nptp.tpr -on SLV.ndx -select "(resname SLV)" """) @@ -596,9 +599,9 @@ quit""",file=fout) if prod == 'nvt': with open('nvtp.mdp', 'w') as f: - print(nvtp_mdp.format(nprod,opt.dt/2000,half_box,half_box,temperature), file=f) + print(nvtp_mdp.format(nprod,opt.dt/2000,rcoulomb,rvdw,temperature), file=f) gmxscript.grompp(f="nvtp.mdp", c="npt.gro", p="system_GMX.top", o="nvtp.tpr",maxwarn=2) - gmxscript.mdrun(deffnm="nvtp", v=True, ntmpi=1) + gmxscript.mdrun(deffnm="nvtp", v=True, ntmpi='{}'.format(opt.nproc)) gmxscript.trjconv(f="nvtp.trr", s="nvtp.tpr", n="mini.ndx", o="finalsystem.pdb") print("Melakukan perhitungan MSD untuk pelarut") os.system("""gmx select -f nvtp.gro -s nvtp.tpr -on SLV.ndx -select "(resname SLV)" """) @@ -609,9 +612,9 @@ quit""",file=fout) os.system("gmx msd -f nvtp.trr -s nvtp.tpr -n {}.ndx -o msd_{}.xvg".format((i[0]+i[1]+i[-1]).upper(),(i[0]+i[1]+i[-1]).upper())) if prod == 'nve': with open('nvep.mdp', 'w') as f: - print(nvep_mdp.format(nprod,opt.dt/2000,half_box,half_box), file=f) + print(nvep_mdp.format(nprod,opt.dt/2000,rcoulomb,rvdw), file=f) gmxscript.grompp(f="nvep.mdp", c="npt.gro", p="system_GMX.top", o="nvep.tpr",maxwarn=2) - gmxscript.mdrun(deffnm="nvep", v=True, ntmpi=1) + gmxscript.mdrun(deffnm="nvep", v=True, ntmpi='{}'.format(opt.nproc)) gmxscript.trjconv(f="nvep.trr", s="nvep.tpr", n="mini.ndx", o="finalsystem.pdb") print("Melakukan perhitungan MSD untuk pelarut") os.system("""gmx select -f nvep.gro -s nvep.tpr -on SLV.ndx -select "(resname SLV)" """) @@ -648,4 +651,4 @@ quit""",file=fout) shutil.rmtree('SistemLarutan', ignore_errors=True) cmmde_solution() # Eksekusi simulasi -make_solution(opt.terlarut,opt.c_terlarut,opt.pelarut,opt.c_pelarut,opt.NumPelarut,opt.persen_terlarut,opt.persen_pelarut,opt.cation,opt.temperature,opt.pressure,opt.packmol,opt.charge_type,opt.lapang,opt.generate_dftbinp,opt.production) +make_solution(opt.terlarut,opt.c_terlarut,opt.pelarut,opt.c_pelarut,opt.NumPelarut,opt.persen_terlarut,opt.persen_pelarut,opt.cation,opt.temperature,opt.pressure,opt.packmol,opt.charge_type,opt.lapang,opt.generate_dftbinp,opt.production,opt.nproc)