#!/usr/local/bin/bash
#CATEGORY: Input creation
#SYNOPSIS: Convert xyz file to optimally-rotated JDFTx ionpos and lattice

if [ -z "$2" ] || [ "$1" == "-h" ] || [ "$1" == "--help" ]; then
	echo '
	Convert an xyz format file (cartesian Angstrom) to a JDFTx geometry specification.
	
	   Usage: xyzToIonposOpt <file>.xyz <pad>
	   
	NOTE: be sure to include command "coords-type cartesian" in the input file!
	
	This will output two files, <file>.ionpos containing Cartesian bohr coordinates
	and <file>.lattice containing the optimum orthogonal cell. The unit cell will be
	larger in each dimension by <pad> bohrs than any atom. The molecule bounding box
	will be centered at the origin, and rotated to minimize the unit cell volume.
	
	This generates the optimum geometry for use with "coulomb-interaction Isolated"
	and "coulomb-truncation-embed 0 0 0". (NOTE: requires octave to be available in path)
	'
	exit 0;
fi

inFile="$1"
ionposFile="${inFile/.xyz/.ionpos}"
latticeFile="${inFile/.xyz/.lattice}"
pad="$2"

cat > $infile.tmp.m <<EOF
if (! exist("fminsearch") )
	pkg load optim #Not loaded by default on many systems
endif

#Read the ion names:
fin = popen("awk 'NR==1 {nAtoms=\$1} NR>2 && NR<=2+nAtoms {print \$1}' $inFile", "r");
nIons = 0;
while(!feof(fin))
	nIons = nIons+1;
	ionnames{nIons} = fgetl(fin);
endwhile
pclose(fin);

#Read the ion coordinates (cartesian Angstrom):
global ionpos;
fin = popen("awk 'NR==1 {nAtoms=\$1} NR>2 && ((NR-1)%(nAtoms+2)>1) {print \$2, \$3, \$4}' $inFile", "r");
ionpos = reshape(fscanf(fin, "%f"), 3,[]);
nConfigs = size(ionpos,2) / nIons
pclose(fin);

ionpos = ionpos / 0.5291772; #Convert to Bohrs

function m = rotAxis(iAxis, theta)
	c = cos(theta); s = sin(theta);
	m = circshift([ c s 0; -s c 0; 0 0 1 ], [iAxis iAxis]);
endfunction

function [out,rot,offs,boxSize] = transformIonpos(in, params)
	rot = rotAxis(3,params(1)) * rotAxis(2,params(2)) * rotAxis(3,params(3)); #ZYZ Euler rotation
	out = rot * in;
	offs = -0.5*(max(out',[],1) + min(out',[],1));
	out += repmat(offs', 1,size(out,2));
	boxSize = max(out',[],1) - min(out',[],1) + $pad;
endfunction

function res = residual(params)
	global ionpos;
	[out,rot,offs,boxSize] = transformIonpos(ionpos, params);
	res = prod(boxSize);
endfunction

initialVolume = residual([0 0 0])
params = fminsearch(@residual, [0 0 0]);
[ionpos,rotationMatrix,offset,boxSize] = transformIonpos(ionpos, params);
optimumVolume = prod(boxSize)
rotationMatrix
offset
boxSize

#Write lattice:
fp = fopen("$latticeFile", "w");
fprintf(fp, "lattice \\\\\n\t%g 0 0 \\\\\n\t0 %g 0 \\\\\n\t0 0 %g\n", boxSize);
fclose(fp);

#Write ionpos file:
for iConfig = 1:nConfigs
	if nConfigs > 1
		fname = sprintf("$ionposFile.%d", iConfig-1);
	else
		fname = "$ionposFile";
	endif
	fp = fopen(fname, "w");
	fprintf(fp, "#Ionic positions in cartesian coordinates:\n");
	for i=1:nIons
		fprintf(fp, "ion %2s %10.5f %10.5f %10.5f 1\n", ionnames{i}, ionpos(:,i+(iConfig-1)*nIons));
	endfor
	fclose(fp);
endfor
EOF

octave -qf $infile.tmp.m
rm $infile.tmp.m
