diff --git a/bin/cmmde.py b/bin/cmmde.py index 35850ae..90c5f64 100755 --- a/bin/cmmde.py +++ b/bin/cmmde.py @@ -60,6 +60,8 @@ parser.add_argument('-l', '--lapang', type=float, default=5, help='Panjang tamba 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.') +parser.add_argument('-nnpt','--nnpt',type=int, default=50000, help='Jumlah step pada saat ekuilibrasi NPT.') +parser.add_argument('-comp','--compress',type=float, default=4.5e-6, help='Nilai kompresibilitas. Default=4.5e-6.') parser.add_argument('-nprod','--nprod',type=int, default=400000, help='Jumlah step pada saat production.') parser.add_argument('-prod', '--production', type=str, default='None', help='Ensembel yang digunakan dalam production run. Pilihan:nve,npt,nvt.') parser.add_argument('-cp', '--c_pelarut', type=int, default=0, help='Muatan bersih molekul pelarut') @@ -373,7 +375,7 @@ if opt.software == 'gromacs': #SBATCH --time=168:0:0 export OMP_NUM_THREADS={} cd $PWD -$GROMACS_COMMAND -mt {} -mp {} -ct {} -pt {} -pp {} -l {} -cat {} -gen {} -Nump {} -prod {} -nprod {} -nequil {} -dt {} -ctype {}""".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),file=fout) +$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) os.system('sbatch run.sh') if opt.job == 'proprep' and opt.software == 'gmx': diff --git a/lib/__pycache__/cmmde_hubbard.cpython-39.pyc b/lib/__pycache__/cmmde_hubbard.cpython-39.pyc index 0cbf362..a6d71b5 100644 Binary files a/lib/__pycache__/cmmde_hubbard.cpython-39.pyc and b/lib/__pycache__/cmmde_hubbard.cpython-39.pyc differ diff --git a/lib/__pycache__/cmmde_orca.cpython-39.pyc b/lib/__pycache__/cmmde_orca.cpython-39.pyc index ffa2e7c..7d0bde6 100644 Binary files a/lib/__pycache__/cmmde_orca.cpython-39.pyc and b/lib/__pycache__/cmmde_orca.cpython-39.pyc differ diff --git a/lib/cmmde_hubbard.py b/lib/cmmde_hubbard.py index c134ab4..49375e8 100644 --- a/lib/cmmde_hubbard.py +++ b/lib/cmmde_hubbard.py @@ -19,6 +19,8 @@ def azimuth(element): 'B' : '2', 'Au' : '3', 'Rh' : '3', + 'Ni' : '3', + 'Co' : '3', 'Ti' : '3', 'C' : '2' } @@ -26,4 +28,4 @@ def azimuth(element): def hubbard(element): U = {'Br':'-0.0573','C':'-0.1492','Ca':'-0.0340','Cl':'-0.0697','F':'-0.1623','H':'-0.1857','I':'-0.0433','K':'-0.0339','Mg':'-0.02','N':'-0.1535','Na':'-0.0454','O':'-0.1575','P':'-0.14','S':'-0.11','Zn':'-0.03','Li': '-0.04548','B':'-0.145'} - return U[element] \ No newline at end of file + return U[element] diff --git a/lib/cmmde_orca.py b/lib/cmmde_orca.py index b517b26..c261592 100644 --- a/lib/cmmde_orca.py +++ b/lib/cmmde_orca.py @@ -5,10 +5,14 @@ def orca(job,method,nproc,geom,charge,mult,scalfreq,temperature,pressure,nroots, method += ' Numfreq' if 'irc' in job: method += ' IRC' - if 'ts' in job: + if 'ts' in job and 'neb' not in job: method += ' OptTS' - if 'NEB-TS' in job: - method += ' NEB-TS' + if 'neb-ts' in job: + method += ' NEB-TS Freq' + if 'neb' in job and 'ts' not in job: + method += ' NEB' + if 'neb-ci' in job: + method += ' NEB-CI' if 'freq' in job and ('XTB2' not in method and 'XTB' not in method and 'XTBFF' not in method): method += ' freq' if disp != 'None': @@ -66,10 +70,10 @@ end""".format(scalfreq,temperature,pressure),file=f) nroots {} tda {} end""".format(nroots,tda),file=f) - if 'NEB' in job: + if 'neb' in job: print("""%NEB -NEB_END_XYZFILE {} -NEB_TS_XYZFILE {} +NEB_END_XYZFILE "{}" +NEB_TS_XYZFILE "{}" END""".format(product,ts),file=f) if 'irc' in job: print("""%irc diff --git a/lib/cmmde_solution.py b/lib/cmmde_solution.py index 201bf11..0c2bf58 100755 --- a/lib/cmmde_solution.py +++ b/lib/cmmde_solution.py @@ -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.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 = 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")