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.

226 lines
11 KiB

2 years ago
import os
from cmmde_hubbard import hubbard, azimuth
def dcdftb(job,method,geom,charge,mult,dispersion,para_path,temp,pressure,ensembel,thermostat,deltat,step,mdprint,a1,a2,a3,b1,b2,b3,c1,c2,c3,restart,traject,velocity,dftbinp,softpot,softtype,softrange,softcenter,metarest,metafreq,metaheight,cvtype,metawidth,pow1,pow2,rcut,fesstart,fesend,fesbin,ag1,ag2,ag3,ag4,solvent,nroots,target,transmult,ocsstr,writetrans,lc,bufrad,delta,opttype,freqtype,econv,dconv):
symb = []
x = []
y = []
z = []
with open(geom, 'r') as f:
Natom = int(next(f).strip())
next(f)
for line in f:
arr = line.split()
symb.append(arr[0])
x.append(arr[1])
y.append(arr[2])
z.append(arr[3])
# Menghitung jumlah atom
elements = list(set(symb))
Numelements = len(elements)
if 'TV' in elements:
elements.remove('TV')
Numelements = len(elements)
with open("cmmd.in",'w') as f:
if method == "DFTB2_gammah":
print("DC=FALSE",file=f)
print("SCC=(LC={} DAMPXH=TRUE DAMPXHZETA=4.5 ECONV={} DCONV={})".format(lc,econv,dconv), file=f)
elif method == "DFTB2":
print("DC=FALSE",file=f)
print("SCC=(LC={} ECONV={} DCONV={})".format(lc, econv, dconv), file=f)
elif method == "DFTB3":
print("DC=FALSE",file=f)
if dispersion != 'D3H5':
print("SCC=(THIRDFULL=TRUE DAMPXH=TRUE DAMPXHZETA=4.0 LC={} ECONV={} DCONV={})".format(lc, econv, dconv),file=f)
else:
print("SCC=(THIRDFULL=TRUE DAMPXH=FALSE LC={} ECONV={} DCONV={})".format(lc, econv, dconv),file=f)
elif method == "DFTB3-diag":
print("DC=FALSE",file=f)
print("SCC=(THIRDDIAG=TRUE DAMPXH=TRUE DAMPXHZETA=4.95 LC={} ECONV={} DCONV={})".format(lc, econv, dconv),file=f)
elif method == "DCDFTB2_gammah":
print("DC=(BUFRAD={} DELTARXYZ={})".format(bufrad,delta),file=f)
print("SCC=(LC={} DAMPXH=TRUE DAMPXHZETA=4.5 ECONV={} DCONV={})".format(lc,econv,dconv), file=f)
elif method == "DCDFTB2":
print("DC=(BUFRAD={} DELTARXYZ={})".format(bufrad,delta),file=f)
print("SCC=(LC={} ECONV={} DCONV={})".foramt(lc, econv, dconv),file=f)
elif method == "DCDFTB3":
print("DC=(BUFRAD={} DELTARXYZ={})".format(bufrad,delta),file=f)
print("SCC=(THIRDFULL=TRUE DAMPXH=TRUE DAMPXHZETA=4.0 LC={} ECONV={} DCONV={})".format(lc, econv, dconv),file=f)
elif method == "DCDFTB3-diag":
print("DC=(BUFRAD={} DELTARXYZ={}))".format(bufrad,delta),file=f)
print("SCC=(THIRDDIAG=TRUE DAMPXH=TRUE DAMPXHZETA=4.95 LC={} ECONV={} DCONV={})".format(lc, econv, dconv),file=f)
# Mapping pelarut implisit
gbsolventtype = { 'water': '1','acetonitril': '2','dimethylsulfoxide':'3','methanol':'4','chloroform':'5','dichloromethane':'6','benzene':'7','toluene':'8','acetone':'9','tetrahydrofuran':'10','diethylether':'11','carbondisulfide':'12','dimethylformamide':'13','hexane':'14','nitromethane':'15','furan':'16','dioxane':'17','ethylacetate':'18','phenol':'19','aniline':'20','benzaldehyde':'21','octanol':'23','dimethylacetamide':'26','ethanol':'25','formamide':'27','methylethylketone':'28'}
if solvent != 'none':
if job == 'td':
if solvent == 'ethanol': gbsolventtype[solvent] = 4
print("SCC=(ALPB=TRUE GBSOLVENTTYPE={} GBPARAMTYPE=7)".format(gbsolventtype[solvent]),file=f)
print("SCC=(SA=TRUE GBSAHBOND=TRUE)",file=f)
else:
print("SCC=(ALPB=TRUE GBSOLVENTTYPE={} GBPARAMTYPE=4)".format(gbsolventtype[solvent]),file=f)
print("SCC=(SA=TRUE GBSAHBOND=TRUE)",file=f)
if job == 'td':
print("TD = (NSTATE={} TARGETSTATE={} MULT={} OSCSTR={} WRITETRANSITION = {})".format(nroots,target,transmult,ocsstr,writetrans),file=f)
if dispersion == "D3":
print("DISPERSION=(DISPTYPE=4)",file=f)
elif dispersion == "D2":
print("DISPERSION=(DISPTYPE=1)",file=f)
elif dispersion == "SK":
print("DISPERSION=(DISPTYPE=2)",file=f)
elif dispersion == "LJ":
print("DISPERSION=(DISPTYPE=3)",file=f)
elif dispersion == "D3":
print("DISPERSION=(DISPTYPE=4)",file=f)
elif dispersion == "D3BJ":
print("DISPERSION=(DISPTYPE=5)",file=f)
elif dispersion == "D3H4":
print("DISPERSION=(DISPTYPE=6)",file=f)
elif dispersion == "DFTulg":
print("DISPERSION=(DISPTYPE=7)",file=f)
elif dispersion == "dDMC":
print("DISPERSION=(DISPTYPE=8)",file=f)
elif dispersion == "D3H5":
print("DISPERSION=(DISPTYPE=9)",file=f)
elif dispersion == "kubilius":
print("DISPERSION=(DISPTYPE=10)",file=f)
elif dispersion == "goursot":
print("DISPERSION=(DISPTYPE=11)",file=f)
elif dispersion == "vdwts":
print("DISPERSION=(DISPTYPE=12)",file=f)
elif dispersion == "manybody":
print("DISPERSION=(DISPTYPE=13)",file=f)
elif dispersion == "D4":
print("DISPERSION=(DISPTYPE=14)",file=f)
if "opt" in job:
opt_type={'bfgs':'1','sd':'2','cg':'3','qm':'4','fire':'5'}
print("OPT=(MAXITER=1000 OPTTYPE={})".format(opt_type[opttype]),file=f)
if "freq" in job:
print("FREQ=(THERMOTEMP={} THERMOPRES={} FREQTYPE={})".format(temp,pressure,freqtype),file=f)
if ensembel == 'NVE':
nvt = 'FALSE'
elif ensembel == 'NVT':
nvt = 'TRUE'
if thermostat == 'nose':
thermo = 3
elif thermostat == 'berendsen':
thermo = 4
elif thermostat == 'andersen':
thermo = 5
if job == "md":
if ensembel == 'NVT':
nvt = 'TRUE'
else:
nvt = 'FALSE'
print('MD=(NVT={} NVTTYPE={} BATHTEMP={} ERRORTEMP=100)'.format(nvt,thermo,temp),file=f)
print('MD=(DELTAT={} NSTEP={} PRINT={})'.format(deltat*1e-15,step,mdprint),file=f)
if softpot == 'true':
print('MD=(SOFT=TRUE SOFTSHAPETYPE={} SOFTRANGE={} SOFTCENTERTYPE={})'.format(softtype,softrange,softcenter))
if restart == 'true':
print('MD=(READVELOCITY=TRUE)',file=f)
if job == "mtd":
if ensembel == 'NVT':
nvt = 'TRUE'
else:
nvt = 'FALSE'
if metarest == 'false':
metarest = 'FALSE'
else:
metarest = 'TRUE'
print('MD=(NVT={} NVTTYPE={} BATHTEMP={} ERRORTEMP=100)'.format(nvt,thermo,temp),file=f)
print('MD=(DELTAT={} NSTEP={} PRINT={})'.format(deltat*1e-15,step,mdprint),file=f)
print('MD=(METADYNAMICS=TRUE METARESTART={} METAPRINTFES=TRUE METAFREQ={} METAHEIGHT={})'.format(metarest,metafreq,metaheight),file=f)
if softpot == 'true':
print('MD=(SOFT=TRUE SOFTSHAPETYPE={} SOFTRANGE={} SOFTCENTERTYPE={})'.format(softtype,softrange,softcenter),file=f)
if restart == 'true':
print('MD=(READVELOCITY=TRUE)',file=f)
with open('metacv.dat','w') as fcv:
atomgroup1 = ag1.split()
atomgroup2 = ag2.split()
nag1 = len(atomgroup1)
nag2 = len(atomgroup2)
if cvtype == 'coordnum':
print("""RATIONALCOORDINATIONNUMBER {} {} {} {} {} {} LABELS AVERAGE {} {} {}
{}
{}
""".format(metawidth,nag1,nag2,pow1,pow2,rcut,fesstart,fesend,fesbin,ag1,ag2),file=fcv)
if cvtype == 'distance':
print("BONDDISTANCE {} {} {} {} {} {}".format(metawidth,ag1,ag2,fesstart,fesend,fesbin),file=fcv)
if cvtype == 'angle':
print("BONDANGLE {} {} {} {} {} {} {}".format(metawidth,ag1,ag2,ag3,fesstart,fesend,fesbin),file=fcv)
if cvtype == 'dihedral':
print("BONDDIHEDRAL {} {} {} {} {} {} {} {}".format(metawidth,ag1,ag2,ag3,ag4,fesstart,fesend,fesbin), file=fcv)
if cvtype == 'distancediff':
print("BONDDISTANCEDIFFERENCE {} {} {} {} {} {} {} {}".format(metawidth,ag1,ag2,ag3,ag4,fesstart,fesend,fesbin),file=fcv)
if cvtype == 'distanceadd':
print("BONDDISTANCEADDITION {} {} {} {} {} {} {} {}".format(metawidth,ag1,ag2,ag3,ag4,fesstart,fesend,fesbin),file=fcv)
if cvtype == 'meandistance':
print("""MEANDISTANCE {} {} {} {} LABELS {} {} {}
{}
{}
""".format(metawidth,nag1,nag2,pow1,fesstart,fesend,fesbin,ag1,ag2),file=fcv)
if cvtype == 'pointplanedistance':
print("ATOMPOINTPLANCEDISTANCE {} {} {} {} {}".format(metawidth,ag1,ag2,ag3,ag4),file=fcv)
print("",file=f)
print("DCDFTB input generated from CMMDE code",file=f)
print("",file=f)
print(Numelements,file=f)
# Pemetaan turunan Hubbard untuk atom-atom tertentu:
for element in elements:
if method == 'DFTB2' or method == 'DCDFTB2' or method == 'DFTB3-diag' or method == 'DCDFTB3-diag' or method == 'DFTB2_gammah' or method == 'DCDFTB2_gammah':
print(element,azimuth(element),file=f)
elif method == 'DFTB3' or method == 'DCDFTB3' or method == 'DCDFTB3-diag':
print(element,azimuth(element),hubbard(element),file=f)
ind = 0
for element2 in elements:
os.system("cp {}/{}-{}.skf .".format(para_path,element,element2))
if ind < len(elements):
print(element+'-'+element2+'.skf',end=' ',file=f)
ind+=1
if ind == len(elements):
print(' ',file=f)
print("",file=f)
print("{} {} {}".format(Natom,charge,mult),file=f)
if restart == 'false':
for symb, x, y, z in zip(symb,x,y,z):
print("{} {} {} {}".format(symb,x,y,z), file=f)
else:
with open(traject, 'r') as ftraj:
lines = ftraj.readlines()
coords = []
for i in range(-Natom,0):
coords.append(lines[i].strip('\n'))
with open(velocity,'r') as fvel:
lines = fvel.readlines()
velocities = []
for i in range(-Natom,0):
velocities.append(lines[i].strip('\n'))
with open('veloc.dat','w') as fveloc:
for i in velocities:
print(i[3:], file=fveloc)
with open(dftbinp, 'r') as inp:
latt_info = []
for line in inp:
if 'TV' in line:
arr = line.strip('\n')
latt_info.append(arr)
for i in coords:
print(i, file=f)
for i in latt_info:
print(i, file=f)
if a1 != 0:
print("""TV {} {} {}
TV {} {} {}
TV {} {} {}""".format(a1,a2,a3,b1,b2,b3,c1,c2,c3),file=f)
print(" ",file=f)