<reponame>aspuru-guzik-group/kraken
#!/usr/bin/env python
import numpy as np
import sys as sys
import argparse as ap
import pathlib as pl
import subprocess as sp
from ded_pdb import pdb
import matplotlib.pyplot as plt


# Define external program directories
# These paths need to be adjusted!
dftd3 = pl.Path("/uufs/chpc.utah.edu/common/home/u1209999/PL_workflow/new_org_use/Pint/dftd3")
dftd4 = pl.Path("/uufs/chpc.utah.edu/common/home/u1209999/PL_workflow/new_org_use/Pint/dftd4")
cwd = pl.Path.cwd()

# Define constants
conv_hartree = 627.509474 # conversion to kcal/mol; taken from https://en.wikipedia.org/wiki/Hartree (20.09.2018)
conv_bohr = 0.52917721067 # conversion from bohr to angstrom; taken from https://en.wikipedia.org/wiki/Bohr_radius (20.09.2018)

def rm_anch(atoms, xyz, surf, anchor, cn=None):
	# Define vdW radii of atoms
	# Written for HCC anchor group (alkynyl)
	vrad = np.column_stack([np.array(['H','C']),np.array([1.20/conv_bohr*1.37,1.77/conv_bohr*1.19])]) # Base values taken from Alvarez 2013, corrections optimized to give clear cuts for Cy and tBu
	l_axis = xyz[anchor[-2]] - xyz[anchor[0]]
	l_norm = np.sqrt(np.sum((l_axis)**2))
	
	surf_del = []
	for i in range(len(surf)):
		# Create cylinder around anchor group and remove atoms
		l0 = abs(np.dot(surf[i] - xyz[anchor[0]],l_axis)/l_norm)
		l1 = abs(np.dot(surf[i] - xyz[anchor[-2]],l_axis)/l_norm)
		b0 = np.sqrt(np.sum((np.cross(surf[i] - xyz[anchor[0]],l_axis))**2))/l_norm
		ref = np.asscalar(vrad[vrad[:,0] == atoms[anchor[-2]],1].astype(float))
		
		if l0 <= l_norm and l1 <= l_norm and b0 <= ref:
			surf_del.append(i)
			continue
		else:
			# Remove all points within certain distance of anchor atoms
			for j in range(len(anchor)-1):
				ij_ref = np.asscalar(vrad[vrad[:,0] == atoms[anchor[j]],1].astype(float)) # get reference distance
				ij_anchor = np.sqrt(np.sum((xyz[anchor[j]] - surf[i])**2))
				ij_attached = np.sqrt(np.sum((xyz[anchor[-1]] - surf[i])**2))
				
				if ij_anchor < ij_ref and ij_attached > ij_anchor:
					surf_del.append(i)
					break
				else:
					continue
	
	for ni in np.sort(surf_del)[::-1]:
		surf = np.delete(surf, ni, 0)
	for ni in np.sort(anchor[:-1])[::-1]:
		atoms = np.delete(atoms, ni, 0)
		xyz = np.delete(xyz, ni, 0)
		if cn != None:
			del cn[ni]

	if cn == None:
		return atoms, xyz, surf
	else:
		return atoms, xyz, surf, cn

def disp(xyz, method, charge=0):
	# Run DFTD3
	if method.casefold() == "d3".casefold():
		d3_out = sp.check_output([dftd3, xyz, "-func", "b3-lyp", "-bjm"])
		d3_out = d3_out.decode('utf-8')
		with open(xyz.split(".")[-2]+"_d3.out","w") as f:
			f.write(d3_out)
		# Read output
		start = "#               XYZ [au]               R0(AA) [Ang.]  CN       C6(AA)     C8(AA)   C10(AA) [au]"
		end = "molecular C6(AA) [au]"
		read = 0
		d3 = []
		for line in d3_out.splitlines():
			if start in line and read == 0:
				read = 1
				continue
			elif end in line and read == 1:
				break
			elif read == 1:
				if line.split() != []:
					d3.append([float(line.split()[7]), float(line.split()[8])])
				else:
					continue
		return d3
	
	# Run DFTD4
	elif method.casefold() == "d4".casefold():
		d4_out = sp.check_output([dftd4, xyz, "-c", str(charge), "--molc6"])
		d4_out = d4_out.decode('utf-8')
		with open(xyz.split(".")[-2]+"_d4.out","w") as f:
			f.write(d4_out)
		# Read output
		start = "#   Z        covCN         q      C6AA      C8AA      α(0)"
		end = " Mol. C6AA /au"
		read = 0
		d4 = []
		for line in d4_out.splitlines():
			if start in line and read == 0:
				read = 1
				continue
			elif end in line and read == 1:
				break
			elif read == 1:
				if line.split() != []:
					d4.append([float(line.split()[5]),float(line.split()[6])])
				else:
					continue
		return d4

def main():
	# Parse arguments
	parser = ap.ArgumentParser()
	parser.add_argument("xyz", help="XYZ file of input structure", type=str)
	parser.add_argument("surf", help="TXT file of points on vdW surface generated by Multiwfn", type=str)
	parser.add_argument("-a", "--anchor", help="Anchor group file with atom numbers to remove. One atom number per line.", type=str, default=None)
	parser.add_argument("-c", "--charge", help="Molecular charge (only for D4).", type=int, default=0)
	parser.add_argument("-d", "--disp", help="Dispersion correction: d3=Grimme D3 Model, d4=Grimme D4 Model", type=str, default="d4")
	parser.add_argument("-n", "--coef", help="Cn coefficients to include: 1=C6, 2=C6+C8, 3=C6+C8+C9(abc)", type=int, default=2)
	args = parser.parse_args()
	
	# Determine file ending
	form = pl.Path(args.xyz).suffix
	if form not in [".xyz"]:
		raise RuntimeError("Unsupported file type.")

	filename = pl.Path(args.xyz).stem

	# Load files
	atoms = np.genfromtxt(pl.Path(args.xyz).resolve(), delimiter=None, usecols = (0), skip_header = 2, dtype=str)
	xyz = np.genfromtxt(pl.Path(args.xyz).resolve(), delimiter=None, usecols = (1, 2, 3), skip_header = 2)/conv_bohr
	surf = np.genfromtxt(pl.Path(args.surf).resolve(), delimiter=None, usecols = (0, 1, 2), skip_header = 1)
	
	# Calculate Cn coefficients
	cn = disp(pl.Path(args.xyz).resolve().name, args.disp, args.charge)
	
	# Remove anchor group
	if args.anchor != None:
		anchor = np.genfromtxt(pl.Path(args.anchor).resolve(), delimiter=None, usecols = (0), skip_header = 0, dtype=int)
		anchor -= 1 # numbering in python starts with 0
		anchor.tolist()
		atoms, xyz, surf, cn = rm_anch(atoms, xyz, surf, anchor, cn)
	natoms = len(xyz)
	npoints = len(surf)
	
	# Save Cn coefficients to file
	with open(cwd / pl.Path(filename+"_ded_c6_" + args.disp + ".out"), 'w') as f:
		for ni in range(natoms):
			f.write(format(str(atoms[ni]), '>3s') + "  " + format(format(cn[ni][0], '.4f'), '>12s') + "  " + format(format(cn[ni][1], '.4f'), '>12s') + "\n")
	
	# Define pint array
	pint = np.zeros(npoints)

	# Perform pint calculations
	for j in range(npoints):
		for i in range(natoms):
			ij_dist = np.sqrt(np.sum((xyz[i] - surf[j])**2))
			pint[j] += (np.sqrt(cn[i][0])/((ij_dist)**3))*np.sqrt(conv_hartree)
			if args.coef > 1:
				pint[j] += (np.sqrt(cn[i][1])/((ij_dist)**4))*np.sqrt(conv_hartree)
				if args.coef > 2:
					pass # include abc dispersion here
	
	pint_ave = np.mean(pint)
	pint_std = np.std(pint)
	pint_max = np.mean(np.sort(pint)[-10:]) # to increase robustness
	pint_min = np.median(np.sort(pint)[:100]) # to increase robustness

	# plt.figure(figsize=(5, 5))
	# hist,bins = np.histogram(pint,bins="auto")#"auto"
	# plt.hist(pint, bins, alpha=0.5, label='pint',color="black")
	# # plt.legend(loc='best')
	# plt.xlabel("pint",fontsize=20)
	# plt.ylabel("N points",fontsize=20)
	# plt.xticks(fontsize=15)
	# plt.yticks(fontsize=15)
	# plt.savefig("pinthist.png")
	
	# Save pint surface data
	with open(cwd / pl.Path(filename+"_ded_surf_" + args.disp + ".out"), 'w') as f:
		for ni in range(npoints):
			f.write(format(format(surf[ni][0]*conv_bohr, '.7f'), '>14s') + "  " + format(format(surf[ni][1]*conv_bohr, '.7f'), '>14s') + "  " + format(format(surf[ni][2]*conv_bohr, '.7f'), '>14s') + "  " + format(format(pint[ni], '.7f'), '>14s') + "\n")
	
	fpdb = pdb(cwd / pl.Path(filename+"_ded_surf_" + args.disp + ".out"), cwd / pl.Path(filename+"_ded_surf_" + args.disp + ".out"), 0, 3)
	fpdb.dump(filename+"_ded_surf_" + args.disp)
	
	# Save pint report
	with open(cwd / pl.Path(filename+"_ded_" + args.disp + ".txt"), 'w') as f:
		f.write('P PARAMETERS')
		f.write('\n')
		f.write('\n')
		f.write('Pint ' + format(format(pint_ave, '.2f'), '>6s'))
		f.write('\n')
		f.write('dP ' + format(format(pint_std, '.2f'), '>6s'))
		f.write('\n')
		f.write('Pmin ' + format(format(pint_min, '.2f'), '>6s'))
		f.write('\n')
		f.write('Pmax ' + format(format(pint_max, '.2f'), '>6s'))
		f.write('\n')

if __name__ == "__main__":
	main()

#!/usr/bin/env python
import numpy as np
import sys as sys
import argparse as ap
import pathlib as pl
import subprocess as sp
from ded_pdb import pdb
import matplotlib.pyplot as plt


# Define external program directories
# These paths need to be adjusted!
dftd3 = pl.Path("/uufs/chpc.utah.edu/common/home/u1209999/PL_workflow/new_org_use/Pint/dftd3")
dftd4 = pl.Path("/uufs/chpc.utah.edu/common/home/u1209999/PL_workflow/new_org_use/Pint/dftd4")
cwd = pl.Path.cwd()

# Define constants
conv_hartree = 627.509474 # conversion to kcal/mol; taken from https://en.wikipedia.org/wiki/Hartree (20.09.2018)
conv_bohr = 0.52917721067 # conversion from bohr to angstrom; taken from https://en.wikipedia.org/wiki/Bohr_radius (20.09.2018)

def rm_anch(atoms, xyz, surf, anchor, cn=None):
	# Define vdW radii of atoms
	# Written for HCC anchor group (alkynyl)
	vrad = np.column_stack([np.array(['H','C']),np.array([1.20/conv_bohr*1.37,1.77/conv_bohr*1.19])]) # Base values taken from Alvarez 2013, corrections optimized to give clear cuts for Cy and tBu
	l_axis = xyz[anchor[-2]] - xyz[anchor[0]]
	l_norm = np.sqrt(np.sum((l_axis)**2))
	
	surf_del = []
	for i in range(len(surf)):
		# Create cylinder around anchor group and remove atoms
		l0 = abs(np.dot(surf[i] - xyz[anchor[0]],l_axis)/l_norm)
		l1 = abs(np.dot(surf[i] - xyz[anchor[-2]],l_axis)/l_norm)
		b0 = np.sqrt(np.sum((np.cross(surf[i] - xyz[anchor[0]],l_axis))**2))/l_norm
		ref = np.asscalar(vrad[vrad[:,0] == atoms[anchor[-2]],1].astype(float))
		
		if l0 <= l_norm and l1 <= l_norm and b0 <= ref:
			surf_del.append(i)
			continue
		else:
			# Remove all points within certain distance of anchor atoms
			for j in range(len(anchor)-1):
				ij_ref = np.asscalar(vrad[vrad[:,0] == atoms[anchor[j]],1].astype(float)) # get reference distance
				ij_anchor = np.sqrt(np.sum((xyz[anchor[j]] - surf[i])**2))
				ij_attached = np.sqrt(np.sum((xyz[anchor[-1]] - surf[i])**2))
				
				if ij_anchor < ij_ref and ij_attached > ij_anchor:
					surf_del.append(i)
					break
				else:
					continue
	
	for ni in np.sort(surf_del)[::-1]:
		surf = np.delete(surf, ni, 0)
	for ni in np.sort(anchor[:-1])[::-1]:
		atoms = np.delete(atoms, ni, 0)
		xyz = np.delete(xyz, ni, 0)
		if cn != None:
			del cn[ni]

	if cn == None:
		return atoms, xyz, surf
	else:
		return atoms, xyz, surf, cn

def disp(xyz, method, charge=0):
	# Run DFTD3
	if method.casefold() == "d3".casefold():
		d3_out = sp.check_output([dftd3, xyz, "-func", "b3-lyp", "-bjm"])
		d3_out = d3_out.decode('utf-8')
		with open(xyz.split(".")[-2]+"_d3.out","w") as f:
			f.write(d3_out)
		# Read output
		start = "#               XYZ [au]               R0(AA) [Ang.]  CN       C6(AA)     C8(AA)   C10(AA) [au]"
		end = "molecular C6(AA) [au]"
		read = 0
		d3 = []
		for line in d3_out.splitlines():
			if start in line and read == 0:
				read = 1
				continue
			elif end in line and read == 1:
				break
			elif read == 1:
				if line.split() != []:
					d3.append([float(line.split()[7]), float(line.split()[8])])
				else:
					continue
		return d3
	
	# Run DFTD4
	elif method.casefold() == "d4".casefold():
		d4_out = sp.check_output([dftd4, xyz, "-c", str(charge), "--molc6"])
		d4_out = d4_out.decode('utf-8')
		with open(xyz.split(".")[-2]+"_d4.out","w") as f:
			f.write(d4_out)
		# Read output
		start = "#   Z        covCN         q      C6AA      C8AA      α(0)"
		end = " Mol. C6AA /au"
		read = 0
		d4 = []
		for line in d4_out.splitlines():
			if start in line and read == 0:
				read = 1
				continue
			elif end in line and read == 1:
				break
			elif read == 1:
				if line.split() != []:
					d4.append([float(line.split()[5]),float(line.split()[6])])
				else:
					continue
		return d4

def main():
	# Parse arguments
	parser = ap.ArgumentParser()
	parser.add_argument("xyz", help="XYZ file of input structure", type=str)
	parser.add_argument("surf", help="TXT file of points on vdW surface generated by Multiwfn", type=str)
	parser.add_argument("-a", "--anchor", help="Anchor group file with atom numbers to remove. One atom number per line.", type=str, default=None)
	parser.add_argument("-c", "--charge", help="Molecular charge (only for D4).", type=int, default=0)
	parser.add_argument("-d", "--disp", help="Dispersion correction: d3=Grimme D3 Model, d4=Grimme D4 Model", type=str, default="d4")
	parser.add_argument("-n", "--coef", help="Cn coefficients to include: 1=C6, 2=C6+C8, 3=C6+C8+C9(abc)", type=int, default=2)
	args = parser.parse_args()
	
	# Determine file ending
	form = pl.Path(args.xyz).suffix
	if form not in [".xyz"]:
		raise RuntimeError("Unsupported file type.")

	filename = pl.Path(args.xyz).stem

	# Load files
	atoms = np.genfromtxt(pl.Path(args.xyz).resolve(), delimiter=None, usecols = (0), skip_header = 2, dtype=str)
	xyz = np.genfromtxt(pl.Path(args.xyz).resolve(), delimiter=None, usecols = (1, 2, 3), skip_header = 2)/conv_bohr
	surf = np.genfromtxt(pl.Path(args.surf).resolve(), delimiter=None, usecols = (0, 1, 2), skip_header = 1)
	
	# Calculate Cn coefficients
	cn = disp(pl.Path(args.xyz).resolve().name, args.disp, args.charge)
	
	# Remove anchor group
	if args.anchor != None:
		anchor = np.genfromtxt(pl.Path(args.anchor).resolve(), delimiter=None, usecols = (0), skip_header = 0, dtype=int)
		anchor -= 1 # numbering in python starts with 0
		anchor.tolist()
		atoms, xyz, surf, cn = rm_anch(atoms, xyz, surf, anchor, cn)
	natoms = len(xyz)
	npoints = len(surf)
	
	# Save Cn coefficients to file
	with open(cwd / pl.Path(filename+"_ded_c6_" + args.disp + ".out"), 'w') as f:
		for ni in range(natoms):
			f.write(format(str(atoms[ni]), '>3s') + "  " + format(format(cn[ni][0], '.4f'), '>12s') + "  " + format(format(cn[ni][1], '.4f'), '>12s') + "\n")
	
	# Define pint array
	pint = np.zeros(npoints)

	# Perform pint calculations
	for j in range(npoints):
		for i in range(natoms):
			ij_dist = np.sqrt(np.sum((xyz[i] - surf[j])**2))
			pint[j] += (np.sqrt(cn[i][0])/((ij_dist)**3))*np.sqrt(conv_hartree)
			if args.coef > 1:
				pint[j] += (np.sqrt(cn[i][1])/((ij_dist)**4))*np.sqrt(conv_hartree)
				if args.coef > 2:
					pass # include abc dispersion here
	
	pint_ave = np.mean(pint)
	pint_std = np.std(pint)
	pint_max = np.mean(np.sort(pint)[-10:]) # to increase robustness
	pint_min = np.median(np.sort(pint)[:100]) # to increase robustness

	# plt.figure(figsize=(5, 5))
	# hist,bins = np.histogram(pint,bins="auto")#"auto"
	# plt.hist(pint, bins, alpha=0.5, label='pint',color="black")
	# # plt.legend(loc='best')
	# plt.xlabel("pint",fontsize=20)
	# plt.ylabel("N points",fontsize=20)
	# plt.xticks(fontsize=15)
	# plt.yticks(fontsize=15)
	# plt.savefig("pinthist.png")
	
	# Save pint surface data
	with open(cwd / pl.Path(filename+"_ded_surf_" + args.disp + ".out"), 'w') as f:
		for ni in range(npoints):
			f.write(format(format(surf[ni][0]*conv_bohr, '.7f'), '>14s') + "  " + format(format(surf[ni][1]*conv_bohr, '.7f'), '>14s') + "  " + format(format(surf[ni][2]*conv_bohr, '.7f'), '>14s') + "  " + format(format(pint[ni], '.7f'), '>14s') + "\n")
	
	fpdb = pdb(cwd / pl.Path(filename+"_ded_surf_" + args.disp + ".out"), cwd / pl.Path(filename+"_ded_surf_" + args.disp + ".out"), 0, 3)
	fpdb.dump(filename+"_ded_surf_" + args.disp)
	
	# Save pint report
	with open(cwd / pl.Path(filename+"_ded_" + args.disp + ".txt"), 'w') as f:
		f.write('P PARAMETERS')
		f.write('\n')
		f.write('\n')
		f.write('Pint ' + format(format(pint_ave, '.2f'), '>6s'))
		f.write('\n')
		f.write('dP ' + format(format(pint_std, '.2f'), '>6s'))
		f.write('\n')
		f.write('Pmin ' + format(format(pint_min, '.2f'), '>6s'))
		f.write('\n')
		f.write('Pmax ' + format(format(pint_max, '.2f'), '>6s'))
		f.write('\n')

if __name__ == "__main__":
	main()
#!/usr/bin/env python # Define external program directories # These paths need to be adjusted! # Define constants # conversion to kcal/mol; taken from https://en.wikipedia.org/wiki/Hartree (20.09.2018) # conversion from bohr to angstrom; taken from https://en.wikipedia.org/wiki/Bohr_radius (20.09.2018) # Define vdW radii of atoms # Written for HCC anchor group (alkynyl) # Base values taken from Alvarez 2013, corrections optimized to give clear cuts for Cy and tBu # Create cylinder around anchor group and remove atoms # Remove all points within certain distance of anchor atoms # get reference distance # Run DFTD3 # Read output # Run DFTD4 # Read output # Parse arguments # Determine file ending # Load files # Calculate Cn coefficients # Remove anchor group # numbering in python starts with 0 # Save Cn coefficients to file # Define pint array # Perform pint calculations # include abc dispersion here # to increase robustness # to increase robustness # plt.figure(figsize=(5, 5)) # hist,bins = np.histogram(pint,bins="auto")#"auto" # plt.hist(pint, bins, alpha=0.5, label='pint',color="black") # # plt.legend(loc='best') # plt.xlabel("pint",fontsize=20) # plt.ylabel("N points",fontsize=20) # plt.xticks(fontsize=15) # plt.yticks(fontsize=15) # plt.savefig("pinthist.png") # Save pint surface data # Save pint report