# 'read_geom' - function for reading coordinates
# input - .xyz or .cosmo file; returns data frame with columns "atom" (element symbol) and "x","y","z" (cartesian coordinates)

read_geom<-function(file) { 
	if (grepl('\\.xyz$',file)) {data<-read.table(file=file, skip=2, col.names=c("atom","x","y","z"))}
	if (grepl('\\.cosmo$',file)) {s<-scan(file=file,what="complex");
		data<-as.data.frame(matrix(s[(grep("!DATE",s)+1):(grep("screening_charge",s)-3)],ncol=9,byrow=TRUE)[,c(8,2,3,4)])
		data[,2:4]<-apply(data[,2:4],2,as.numeric)
		names(data)<-c("atom","x","y","z")}
return(data)}

# function 'len()' - returns lengths of vectors. 
# input -  vector coordinates c(x,y,z) or coordinate matrix with x, y, z as columns

len<-function(xyz){if(is.null(dim(xyz))){xyz<-matrix(xyz,ncol=3)}; return(sqrt(rowSums(xyz^2)))}

# 'min_function' - minimization function for cone angle calculation. Based on Eq. 32 from Bilbrey et al (doi.org/10.1002/jcc.23217)
# input - vector of cone axis ('param'); returns 1/2 of cone angle. 
# NB! Apex is required to be at the origin (0,0,0)

min_function<-function(param){
	rows=rowSums(coord)!=0    # if apex is included in the coordinate matrix, exclude it
	m=apply(X=as.matrix(coord[rows,]),MARGIN=2, function(x) x/len(coord[rows,])) # unit vectors directed from from apex to atom centers
	n=matrix(param,nrow=3)/len(param) # unit vector of cone axis
	beta=asin(data$r[rows]/len(coord[rows,])) # angles between vectors m and vectors tangent to vdW surfaces
	# If apex lies within vdW shpere another atom (e.g. P), the beta-angle such atom does not exist and NA is produced 
	# Replacing it with the highest of the remaining beta-angles solves the problem:
	beta[is.na(beta)]=max(beta,na.rm=TRUE) 
	angles=beta+acos(m%*%n) # half-cone angles
return(max(angles))}

# vdW radii (values used in Bilbrey et al, doi.org/10.1002/jcc.23217):
radii<-c(P=1.8, H=1.2, C=1.7, N=1.55, O=1.52, F=1.47, Cl=1.75, S=1.8)

coordinate_file<-"" # insert path to your file

data<-read_geom(file=coordinate_file)  # read data from file
data$r<-radii[data$atom] # add radii to the input data
coord<-as.matrix(data[,c("x","y","z")]) # coordinate matrix
D<-as.matrix(dist(x=coord, diag=TRUE, upper=TRUE))  # distance matrix for atoms
prot<-which(D[data$atom=="P",]<1.5 & D[data$atom=="P",]>0 & data$atom=="H") # number of the hydrogen atom located within 1.5 A from phosphorus
coord<-sweep(coord,2,coord[prot,]) # shift apex (H) to the origin

	# Optional: remove proton from P and define a virtual apex on continuation on P-H bond
	dist_P_apex<-2.28 # choose apex distance from P atom (in Angstroms)
	apex=coord[data$atom=="P",]-coord[data$atom=="P",]/len(coord[data$atom=="P",])*dist_P_apex  # coordinates of virtual apex
	coord<-sweep(coord,2,apex)[-prot,] # shift apex to origin and remove proton
	data<-data[-prot,] # remove proton from input data

guess<-apply(coord,2,median) 	# initial guess of the cone axis
O<-optim(par=as.numeric(guess)/len(guess), fn=min_function) 	# optimize cone axis for minimum cone angle

O$value*2  # exact cone angle in radians

O$value*360/pi  # exact cone angle in degrees


# visualization:
library(rgl)
col<-c(C="grey",H="white",O="red",N="blue",P="orange",F="blue",Cl="green") # atom colours
coord_ranges<-apply(coord,2,max)-apply(coord,2,min) # coordinate ranges (for setting the aspect ratio)
plot3d(coord, col=col[data$atom], type="s", decorate=FALSE, radius=data$r, aspect=coord_ranges/max(coord_ranges)) # show molecule
spheres3d(c(0,0,0), col="magenta",  radius=0.5)  # show apex