@ -22,20 +22,23 @@ parser.add_argument('-cat','--cation',type=str, default='none',help='Kation yang
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 = ' /home/adit/opt/oldpackmol/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 = steep ; Algorithm ( steep = steepest descent minimization )
emtol = 100.0 ; Stop minimization when the maximum force < 1000.0 kJ / mol / nm
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
@ -50,7 +53,7 @@ pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions
"""
nve_mdp = """
define = - DFLEXIBLE
integrator = md ; leap - frog integrator
integrator = md
nsteps = { } ; 1 * 500000 = 50 ps
dt = 0.001 ; 0.5 fs
nstxout = 100 ; save coordinates every 1.0 ps
@ -81,21 +84,21 @@ pcoupl = no ; no pressure coupling in NVE
; Periodic boundary conditions
pbc = xyz ; 3 - D PBC
; Velocity generation
gen_vel = no ; Velocity generation is off
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 = { } ; 1 * 500000 = 50 ps
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 = no ; first dynamics run
continuation = yes ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h - bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
@ -114,16 +117,14 @@ 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
tau_t = 0.0 1 ; 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 = yes ; assign velocities from Maxwell distribution
gen_temp = { } ; temperature for Maxwell distribution
gen_seed = - 1 ; generate a random seed
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.
@ -131,8 +132,8 @@ npt_mdp = """
title = NPT equilibration
; Run parameters
integrator = md ; leap - frog integrator
nsteps = { } ; 2 * 50000 = 100 ps
dt = 0.002 ; 2 fs
nsteps = { }
dt = 0.001 ; 1 fs
; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
@ -143,7 +144,7 @@ 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
lincs_order = 6 ; also related to accuracy
; Neighborsearching
cutoff - scheme = Verlet
ns_type = grid ; search neighboring grid cells
@ -162,9 +163,9 @@ 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.0 ; time constant , in ps
tau_p = 2 ; time constant , in ps
ref_p = { } ; reference pressure , in bar
compressibility = 4.5e-5 ; isothermal compressibility of water , bar ^ - 1
compressibility = { } ; isothermal compressibility of water , bar ^ - 1
refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3 - D PBC
@ -204,10 +205,10 @@ 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.0 ; time constant , in ps
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
compressibility = 4.5e-5 ; isothermal compressibility of water , bar ^ - 1
refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3 - D PBC
@ -556,16 +557,18 @@ quit""",file=fout)
# 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 , temperature ) , file = f )
print ( nvt_mdp . format ( nequil , half_box , half_box , temperature ) , file = f )
with open ( ' npt.mdp ' , ' w ' ) as f :
print ( npt_mdp . format ( nequil , half_box , half_box , temperature , pressure ) , file = f )
print ( npt_mdp . format ( nnpt , half_box , half_box , temperature , pressure , compress ) , file = f )
gmxscript . grompp ( f = " minim.mdp " , c = " system_box.pdb " , o = " mini.tpr " , p = " system_GMX.top " , maxwarn = 2 )
gmxscript . mdrun ( deffnm = " mini " , ntmpi = 1 )
@ -579,7 +582,7 @@ quit""",file=fout)
if prod == ' npt ' :
with open ( ' nptp.mdp ' , ' w ' ) as f :
print ( nptp_mdp . format ( nprod , opt . dt / 2000 , half_box , half_box , temperature , pressure ) , file = 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 " )
