A series of Python3 script to lower the barrier of computing and simulating molecular and material systems.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 

651 lines
29 KiB

#!/usr/bin/env python3
import subprocess
import os
import shutil
import gmxscript
from scipy.constants import N_A
import argparse
import sys
import time
from cmmde_mass import mass
parser = argparse.ArgumentParser(description='SPE Solution Program')
parser.add_argument('-mt', '--terlarut', type=str, help='Nama file molekul terlarut dalam format .xyz (ditulis tanpa ekstensi).')
parser.add_argument('-ctype','--charge_type',type=str, default='gas', help='Tipe muatan yang digunakan untuk parameterisasi muatan. Pilihan: gas, bcc, dftb')
parser.add_argument('-ct', '--c_terlarut', type=str, help='Muatan bersih molekul terlarut. Jika terlarut lebih dari satu buah, pisahkan dengan koma')
parser.add_argument('-pt','--persen_terlarut',type=str, help='Fraksi massa terlarut dalam persen. Jika pelarut lebih dari satu buah, pisahkan menggunakan koma.')
parser.add_argument('-pp','--persen_pelarut',type=float, help='Fraksi massa pelarut dalam persen.')
parser.add_argument('-mp', '--pelarut', type=str, help='Nama file molekul pelarut dalam format .xyz (ditulis tanpa ekstensi).')
parser.add_argument('-cp', '--c_pelarut', type=int, help='Muatan bersih molekul pelarut',default=0)
parser.add_argument('-Nump','--NumPelarut', type=int, default=100, help='Jumlah molekul pelarut maksimum dalam sistem larutan. Default = 100.')
parser.add_argument('-cat','--cation',type=str, default='none',help='Kation yang digunakan dalam sistem elektrolit baterai.')
parser.add_argument('-p', '--pressure', default=1.0, type=float, help='Tekanan dalam satuan bar')
parser.add_argument('-t', '--temperature', default=298.15, type=float, help='Suhu yang digunakan dalam simulasi.')
parser.add_argument('-nequil','--nequil',type=int, default=50000, help='Jumlah step pada saat ekuilibrasi.')
parser.add_argument('-nnpt','--nnpt',type=int, default=50000, help='Jumlah step pada saat ekuilibrasi NPT.')
parser.add_argument('-nprod','--nprod',type=int, default=400000, help='Jumlah step pada saat production.')
parser.add_argument('-dt','--dt',type=float, default=1, help='Step simulasi dinamika molekul yang dilakukan.')
parser.add_argument('-comp','--compress',type=float, default=4.5e-6, help='Nilai kompresibilitas. Default=4.5e-6.')
parser.add_argument('--packmol', type=str, default='packmol')
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.')
opt=parser.parse_args(sys.argv[1:])
# Input file minimisasi energi dalam GROMACS
minim_mdp = """
integrator = cg ; Algorithm (steep or cg)
emtol = 1 ; Stop minimization when the maximum force < 1 kJ/mol/nm
emstep = 0.01 ; Minimization step size
nsteps = 50000 ; Maximum number of (minimization) steps to perform
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 5 ; Frequency to update the neighbor list and long range forces
ns_type = grid ; Method to determine neighbor list (simple, grid)
coulombtype = PME ; Treatment of long range electrostatic interactions
rcoulomb = {} ; Short-range electrostatic cut-off
rvdw = {} ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions
;define = -DFLEXIBLE
"""
nve_mdp = """
define = -DFLEXIBLE
integrator = md
nsteps = {} ; 1 * 500000 = 50 ps
dt = 0.001 ; 0.5 fs
nstxout = 100 ; save coordinates every 1.0 ps
nstenergy = 100 ; save energies every 1.0 ps
nstlog = 100 ; update log file every 1.0 ps
; Bond parameters
continuation = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet
rcoulomb = {} ; short-range electrostatic cutoff (in nm)
rvdw = {} ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is off
tcoupl = no ; no temperature coupling in NVE
; Pressure coupling is off
pcoupl = no ; no pressure coupling in NVE
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = yes ; Velocity generation is off
"""
# Input file simulasi menggunakan ensembel kanonik (NVT) dalam GROMACS
nvt_mdp="""
define = -DFLEXIBLE
integrator = md ; leap-frog integrator
nsteps = {}
dt = 0.001 ; 1 fs
nstxout = 100 ; save coordinates every 1.0 ps
;nstenergy = 100 ; save energies every 1.0 ps
;nstlog = 100 ; update log file every 1.0 ps
; Bond parameters
continuation = yes ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet
rcoulomb = {} ; short-range electrostatic cutoff (in nm)
rvdw = {} ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = System ; two coupling groups - more accurate
tau_t = 0.01 ; time constant, in ps
ref_t = {} ; reference temperature, one for each group, in K
; Pressure coupling is off
pcoupl = no ; no pressure coupling in NVT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = no ; assign velocities from Maxwell distribution
"""
# Input file simulasi isoterm-isobarik (NPT) dalam GROMACS. Harapannya, setelah dilakukan NPT, akan diperoleh sistem dengan massa jenis yang optimal.
npt_mdp = """
title = NPT equilibration
; Run parameters
integrator = md ; leap-frog integrator
nsteps = {}
dt = 0.001 ; 1 fs
; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
; Bond parameters
continuation = yes ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 6 ; also related to accuracy
; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = {} ; short-range electrostatic cutoff (in nm)
rvdw = {} ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = System ; two coupling groups - more accurate
tau_t = 0.1 ; time constant, in ps
ref_t = {} ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = Berendsen ; Pressure coupling on in NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2 ; time constant, in ps
ref_p = {} ; reference pressure, in bar
compressibility = {} ; isothermal compressibility of water, bar^-1
refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off
"""
nptp_mdp = """
title = NPT Production
; Run parameters
integrator = md ; leap-frog integrator
nsteps = {} ;
dt = {} ; in ps
; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
; Bond parameters
continuation = yes ; if no => first dynamics run
; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = {} ; short-range electrostatic cutoff (in nm)
rvdw = {} ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is off
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = System ; two coupling groups - more accurate
tau_t = 0.1 ; time constant, in ps
ref_t = {} ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = Berendsen ; Pressure coupling on in NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 10 ; time constant, in ps
compressibility = {} ; isothermal compressibility of water, bar^-1
ref_p = {} ; reference pressure, in bar
refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off
"""
nvtp_mdp = """
title = NVT Production
; Run parameters
integrator = md ; leap-frog integrator
nsteps = {} ;
dt = {} ; in ps
; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
; Bond parameters
continuation = yes ; if no => first dynamics run
; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = {} ; short-range electrostatic cutoff (in nm)
rvdw = {} ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is off
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = System ; two coupling groups - more accurate
tau_t = 0.1 ; time constant, in ps
ref_t = {} ; reference temperature, one for each group, in K
; Pressure coupling is off in NVT simulation
pcoupl = no ; Pressure coupling on in NPT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off
"""
nvep_mdp = """
title = NVE Production
; Run parameters
integrator = md ; leap-frog integrator
nsteps = {} ;
dt = {} ; in ps
; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
; Bond parameters
continuation = yes ; if no => first dynamics run
; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = {} ; short-range electrostatic cutoff (in nm)
rvdw = {} ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is off
tcoupl = no ; modified Berendsen thermostat
; Pressure coupling is off in NVE simulation
pcoupl = no ; Pressure coupling on in NPT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off
"""
def chargeparm(mol2file,dftbdat):
header = []
footer = []
serial = []
symbol = []
x = []
y = []
z = []
type = []
resid = []
resname = []
charge = []
with open(mol2file, 'r') as f:
for i in range(1,9):
header.append(next(f).strip())
for line in range(1,int(header[2].split()[0])+1):
arr = next(f).split()
serial.append(arr[0])
symbol.append(arr[1])
x.append(arr[2])
y.append(arr[3])
z.append(arr[4])
type.append(arr[5])
resid.append(arr[6])
resname.append(arr[7])
for line in f:
footer.append(line.strip())
with open(dftbdat,'r') as f:
for line in f:
if 'MULLIKEN' in line:
next(f)
for i in range(1,int(header[2].split()[0])+1):
arr = next(f).split()
charge.append(arr[2])
fname = mol2file.split(".")[0]
with open("{}_mod.mol2".format(fname),'w') as f:
for line in header[0:4]:
print(line,file=f)
print("DFTBCharge",file=f)
print("""
@<TRIPOS>ATOM""",file=f)
for serial, symbol, x, y, z, type, resid, resname, charge in zip(serial, symbol, x, y, z, type, resid, resname, charge):
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):
if temperature:
print (f'Suhu: {temperature:12.4f} K')
# Penyiapan box molekul
terlarut = terlarut.split(",")
c_terlarut = c_terlarut.split(",")
c_terlarut = dict(zip(terlarut,c_terlarut))
for i in terlarut:
os.chdir('{}'.format(i))
if charge_type == 'dftb':
os.system("obabel cmmd.xyz -O {}.pdb".format(i))
else:
os.system("rm geom.xyz")
os.system("rm cmmd_trj.xyz")
os.system("obabel *.xyz -O {}.pdb".format(i))
with open('{}.pdb'.format(i),'r') as f:
fdata = f.read()
fdata = fdata.replace('UNL','{}'.format((i[0]+i[1]+i[-1]).upper()))
fdata = fdata.replace('CONECT', '#CONECT')
with open('{}.pdb'.format(i),'w') as f:
f.write(fdata)
if opt.charge_type == 'dftb':
charge_type='gas'
os.system("antechamber -j 5 -at gaff -dr no -i {}.pdb -fi pdb -o {}.mol2 -fo mol2 -c {} -s 2 -nc {}".format(i,i,charge_type,float(c_terlarut[i])))
chargeparm("{}.mol2".format(i), "dftb.dat")
else:
os.system("antechamber -j 5 -at gaff -dr no -i {}.pdb -fi pdb -o {}.mol2 -fo mol2 -c {} -s 2 -nc {}".format(i,i,charge_type,float(c_terlarut[i])))
os.system("cp {}.mol2 {}_mod.mol2".format(i,i))
os.system("parmchk2 -i {}_mod.mol2 -f mol2 -o {}.frcmod".format(i,i))
with open("tleap.in", 'w') as fout:
print("""
source leaprc.gaff
{} = loadmol2 {}_mod.mol2
loadamberparams {}.frcmod
check {}
saveoff {} {}.lib
saveamberparm {} {}.parmtop {}.inpcrd
savepdb {} {}_gaff.pdb
quit
""".format((i[0]+i[1]+i[-1]).upper(), i, i, (i[0]+i[1]+i[-1]).upper(),(i[0]+i[1]+i[-1]).upper(),(i[0]+i[1]+i[-1]).upper(),(i[0]+i[1]+i[-1]).upper(),(i[0]+i[1]+i[-1]).upper(),(i[0]+i[1]+i[-1]).upper(),(i[0]+i[1]+i[-1]).upper(),(i[0]+i[1]+i[-1]).upper()), file=fout)
os.system("tleap -f tleap.in")
os.chdir('../')
os.chdir('{}'.format(pelarut))
# Penyiapan molekul pelarut menggunakan GAFF
# Memanggil antechamber untuk membuatkan file mol2 yang mengandung informasi muatan bersih atom
if charge_type == 'dftb':
os.system("obabel cmmd.xyz -O {}.pdb".format(pelarut))
else:
os.system("rm geom.xyz")
os.system("rm cmmd_trj.xyz")
os.system("obabel *.xyz -O {}.pdb".format(pelarut))
with open('{}.pdb'.format(pelarut),'r') as f:
fdata = f.read()
fdata = fdata.replace('UNL','SLV')
fdata = fdata.replace('HOH','SLV')
fdata = fdata.replace('CONECT', '#CONECT')
with open('{}.pdb'.format(pelarut),'w') as f:
f.write(fdata)
if opt.charge_type == 'dftb':
charge_type='gas'
os.system("antechamber -j 5 -at gaff -dr no -i {}.pdb -fi pdb -o {}.mol2 -fo mol2 -c {} -s 2 -nc {}".format(pelarut,pelarut,charge_type,c_pelarut))
chargeparm("{}.mol2".format(pelarut), "dftb.dat")
# Mengidentifikasi parameter yang tidak terdefinisikan dalam GAFF
else:
os.system("antechamber -j 5 -at gaff -dr no -i {}.pdb -fi pdb -o {}.mol2 -fo mol2 -c {} -s 2 -nc {}".format(pelarut,pelarut,charge_type,c_pelarut))
os.system("cp {}.mol2 {}_mod.mol2".format(pelarut,pelarut))
os.system("parmchk2 -i {}_mod.mol2 -f mol2 -o {}.frcmod".format(pelarut,pelarut))
with open("tleap.in", 'w') as fout:
print("""
source leaprc.gaff
SLV = loadmol2 {}_mod.mol2
loadamberparams {}.frcmod
check SLV
saveoff SLV {}.lib
saveamberparm SLV {}.parmtop {}.inpcrd
savepdb SLV {}_gaff.pdb
quit
""".format(pelarut,pelarut,pelarut,pelarut,pelarut,pelarut), file=fout)
os.system("tleap -f tleap.in")
if cation != 'none':
os.system("mkdir ../{}".format(cation))
os.chdir('../{}'.format(cation))
with open('{}.mol2'.format(cation),'w') as f:
print("""@<TRIPOS>MOLECULE
ION
1 0 0 0 0
SMALL
GASTEIGER
@<TRIPOS>ATOM
1 {}+ 0.0000 0.0000 0.0000 {}+ 1 ion 1.0000
@<TRIPOS>BOND
""".format(cation,cation),file=f)
# Penyiapan box
os.chdir('../')
with open ("{}/cmmd.xyz".format(pelarut), 'r') as f:
elements = []
Mr_pelarut = 0
next(f)
next(f)
for line in f:
arr = line.split()
elements.append(arr[0])
for ele in elements:
Mr_pelarut += float(mass(ele))
Mr_terlarut = []
for i in terlarut:
ele_terlarut = []
Mr_t = 0
with open("{}/cmmd.xyz".format(i), 'r') as f:
next(f)
next(f)
for line in f:
arr = line.split()
ele_terlarut.append(arr[0])
for ele in ele_terlarut:
Mr_t += float(mass(ele))
Mr_terlarut.append(Mr_t)
Mr_terlarut = dict(zip(terlarut,Mr_terlarut))
m_pelarut = N_pelarut*(Mr_pelarut/N_A) # dalam g
m_larutan = m_pelarut/persen_pelarut*100
persen_terlarut = persen_terlarut.split(",")
persen_terlarut = dict(zip(terlarut, persen_terlarut))
N_terlarut = []
for i in terlarut:
N = round((float(persen_terlarut[i])/100)*m_larutan/Mr_terlarut[i]*N_A)
N_terlarut.append(N)
N_terlarut = dict(zip(terlarut,N_terlarut))
V_box = m_larutan # Dalam satuan cm3
# Panjang rusuk kubus dalam satuan cm:
box_size = V_box**(1/3)
# Panjang rusuk kubus dalam satuan Angstroem:
cmToAng = 1e8
box_size = box_size * cmToAng + lapang
# Hentikan perhitungan jika ukuran kubus terlalu kecil
half_box = min(1.0, box_size/20.0-0.1)
print(half_box)
def cmmde_solution():
cwd = os.path.abspath(os.curdir)
print('Membuat sistem larutan menggunakan PACKMOL')
os.makedirs('SistemLarutan', exist_ok=True)
os.chdir('SistemLarutan')
with open('packmol.inp', 'w') as fout:
print("""tolerance 2.5
filetype pdb
output system_init.pdb
structure ../{}/{}_gaff.pdb
number {}
inside cube -{} -{} -{} {}
resnumbers 3
end structure""".format(pelarut, pelarut, N_pelarut, box_size/2.0, box_size/2.0, box_size/2.0, box_size), file=fout)
for i in terlarut:
print("""structure ../{}/{}_gaff.pdb
number {}
inside cube -{} -{} -{} {}
resnumbers 3
end structure""".format(i,(i[0]+i[1]+i[-1]).upper(),N_terlarut[i],box_size/2.0, box_size/2.0, box_size/2.0, box_size), file=fout)
print('='*80)
with open ('tleap.in', 'w') as fout:
print("""source leaprc.gaff
loadamberparams /Users/adit/opt/amber22/dat/leap/parm/frcmod.ionslm_iod_opc3""", file=fout)
for i in terlarut:
print("""loadoff ../{}/{}.lib
loadamberparams ../{}/{}.frcmod""".format(i,(i[0]+i[1]+i[-1]).upper(),i,i),file=fout)
print("""loadoff ../{}/{}.lib
loadamberparams ../{}/{}.frcmod
SYSTEM = loadpdb system_init.pdb""".format(pelarut,pelarut,pelarut,pelarut),file=fout)
if cation != 'none':
print("""ION = loadmol2 ../{}/{}.mol2
addions SYSTEM ION 0""".format(cation,cation),file=fout)
print("""list
saveamberparm SYSTEM system.prmtop system.inpcrd
savepdb SYSTEM system.pdb
quit""",file=fout)
with open('packmol.inp', 'r') as f:
print(f.read())
print('='*80)
os.system("{} < packmol.inp".format(packmol))
os.system("tleap -f tleap.in")
os.system("acpype -p system.prmtop -x system.inpcrd")
gmxscript.editconf(
f = "system.pdb",
box = [box_size/10.0, box_size/10.0, box_size/10.0],
bt = "cubic",
o = "system_box.pdb"
)
# Menuliskan semua input file yang dibutuhkan, meliputi simulasi NVE, NVT, dan NPT
nequil = opt.nequil
nnpt = opt.nnpt
nprod = opt.nprod
compress = opt.compress
with open('minim.mdp', 'w') as f:
print(minim_mdp.format(half_box, half_box), file=f)
with open('nve.mdp', 'w') as f:
print(nve_mdp.format(nequil,half_box,half_box), file=f)
with open('nvt.mdp', 'w') as f:
print(nvt_mdp.format(nequil,half_box, half_box,temperature), file=f)
with open('npt.mdp', 'w') as f:
print(npt_mdp.format(nnpt,half_box, half_box,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.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.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.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)
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)
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.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)" """)
os.system("gmx msd -f nptp.trr -s nptp.tpr -n SLV.ndx -o msd_SLV.xvg")
print("Melakukan perhitungan MSD untuk terlarut, masukkan masing-masing terlarut")
for i in terlarut:
os.system("""gmx select -f nptp.gro -s nptp.tpr -on {}.ndx -select "(resname {})" """.format((i[0]+i[1]+i[-1]).upper(),(i[0]+i[1]+i[-1]).upper()))
os.system("gmx msd -f nptp.trr -s nptp.tpr -n {}.ndx -o msd_{}.xvg".format((i[0]+i[1]+i[-1]).upper(),(i[0]+i[1]+i[-1]).upper()))
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)
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.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)" """)
os.system("gmx msd -f nvtp.trr -s nvtp.tpr -n SLV.ndx -o msd_SLV.xvg")
print("Melakukan perhitungan MSD untuk terlarut, masukkan masing-masing terlarut")
for i in terlarut:
os.system("""gmx select -f nvtp.gro -s nvtp.tpr -on {}.ndx -select "(resname {})" """.format((i[0]+i[1]+i[-1]).upper(),(i[0]+i[1]+i[-1]).upper()))
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)
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.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)" """)
os.system("gmx msd -f nvep.trr -s nvep.tpr -n SLV.ndx -o msd_SLV.xvg")
print("Melakukan perhitungan MSD untuk terlarut, masukkan masing-masing terlarut")
for i in terlarut:
os.system("""gmx select -f nvep.gro -s nvep.tpr -on {}.ndx -select "(resname {})" """.format((i[0]+i[1]+i[-1]).upper(),(i[0]+i[1]+i[-1]).upper()))
os.system("gmx msd -f nvep.trr -s nvep.tpr -n {}.ndx -o msd_{}.xvg".format((i[0]+i[1]+i[-1]).upper(),(i[0]+i[1]+i[-1]).upper()))
if generate_dftbinp != 'false':
os.system("obabel finalsystem.pdb -O finalsystem.xyz")
with open('finalsystem.xyz', 'r') as f:
fdata = f.read()
fdata = fdata.replace('*','{}'.format(opt.cation))
with open('finalsystem.pdb', 'r') as f:
tv1 = 0
tv2 = 0
tv3 = 0
for line in f:
if 'CRYST' in line:
arr = line.split()
tv1 += float(arr[1])
tv2 += float(arr[2])
tv3 += float(arr[3])
break
with open('finalsystem.xyz'.format(terlarut),'w') as fout:
fout.write(fdata)
print("TV {} 0.00 0.00".format(tv1),file=fout)
print("TV 0.00 {} 0.00".format(tv2), file=fout)
print("TV 0.00 0.00 {}".format(tv3), file=fout)
os.chdir(cwd)
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)