# A script to calculate RESP charge based on ORCA and Multiwfn
# Written by Tian Lu (sobereva@sina.com), last update: 2022-Mar-8
# Examples:
# Calculating neutral singlet molecule in water: RESP_ORCA.sh H2O.xyz
# Calculating anionic singlet molecule in water: RESP_ORCA.sh ani.pdb -1 1
# Calculating neutral triplet molecule in ethanol: RESP_ORCA.sh nico.xyz 0 3 ethanol
# Calculating neutral singlet molecule in ETHYL ETHANOATE: RESP_ORCA.sh nico.xyz 0 1 "ETHYL ETHANOATE"
# Calculating cation doublet molecule in vacuum: RESP_ORCA.sh maki.xyz 0 2 gas

#!/bin/bash

#Set actual paths of ORCA and orca_2mkl utility here
ORCA="/root/orca_6_0_0_shared_openmpi416_avx2/orca"
orca_2mkl="/root/orca_6_0_0_shared_openmpi416_avx2/orca_2mkl"

#Set number of CPU cores used in calculation here
nprocs=12
maxcore=1000

keyword_opt="! B97-3c opt"
keyword_SP="! B3LYP/G D3 def2-TZVP def2/J RIJCOSX"

export inname=$1
filename=${inname%.*}
suffix=${inname##*.}

### Parse arguments
if [ $2 ];then
	echo "Net charge = $2"
	chg=$2
else
	echo "Net charge was not defined. Default to 0"
	chg=0
fi

if [ $3 ];then
	echo "Spin multiplicity = $3"
	multi=$3
else
	echo "Spin multiplicity was not defined. Default to 1"
	multi=1
fi

if [ "$4" ];then
	if [[ $4 == gas ]] ; then
  		echo Calculations will be done in vacuum
  		solvent=""
	else
		echo Solvent is $4
  	solvent=$4
	fi
else
	echo "Solvent name was not defined. Default to water"
  solvent="Water"
fi

### Convert current file to tmp.xyz
Multiwfn_noGUI $1 > /dev/null << EOF
100
2
2
tmp.xyz
0
q
EOF

### Create input file for optimization (opt.inp)
cat << EOF > opt.inp
$keyword_opt
%maxcore $maxcore
%pal nprocs $nprocs end
%geom Convergence loose end
EOF
if [[ $solvent != "" ]] ; then
cat << EOF >> opt.inp
%cpcm
smd true
SMDsolvent "$solvent"
end
EOF
fi
echo "* xyz $chg $multi" >> opt.inp
awk '{if (NR>2) print }' tmp.xyz >> opt.inp
echo "*" >> opt.inp
rm tmp.xyz

### Run optimization
echo Running optimization task via ORCA...
$ORCA opt.inp > opt.out

if grep -Fq "ORCA TERMINATED NORMALLY" opt.out
then
	echo Done!
	mv opt.xyz tmp.xyz
	rm opt.* opt_*
else
	echo The optimization task has failed! Please check content of opt.out to find reason
	echo The script is terminated
	mv opt.out tmp.out
	rm opt.* opt_*
	mv tmp.out opt.out
	exit 1
fi

### Create input file for single point (SP.inp)
cat << EOF > SP.inp
$keyword_SP
%maxcore $maxcore
%pal nprocs $nprocs end
EOF
if [[ $solvent != "" ]] ; then
cat << EOF >> SP.inp
%cpcm
smd true
SMDsolvent "$solvent"
end
EOF
fi
echo "* xyz $chg $multi" >> SP.inp
awk '{if (NR>2) print }' tmp.xyz >> SP.inp
echo "*" >> SP.inp
rm tmp.xyz

### Run single point
echo Running single point task via ORCA...
$ORCA SP.inp > SP.out

if grep -Fq "ORCA TERMINATED NORMALLY" SP.out
then
	echo Done!
else
	echo The single point task has failed! Please check content of SP.out to find reason
	echo The script is terminated
	mv SP.out tmp.out
	rm SP.* SP_*
	mv tmp.out SP.out
	exit 1
fi

### Convert to .molden file
echo Running orca_2mkl...
$orca_2mkl SP -molden > /dev/null

#Field containing number of valence electrons of def2 pseudopotential basis set
cat << EOF > Nval.txt
[Nval]
Rb  9
Sr 10
Y  11
Zr 12
Nb 13
Mo 14
Tc 15
Ru 16
Rh 17
Pd 18
Ag 19
Cd 20
In 21
Sn 22
Sb 23
Te 24
I  25
Xe 26
Cs  9
Ba 10
La 11
Ce 30
Pr 31
Nd 32
Pm 33
Sm 34
Eu 35
Gd 36
Tb 37
Dy 38
Ho 39
Er 40
Tm 41
Yb 42
Lu 43
Hf 12
Ta 13
W  14
Re 15
Os 16
Ir 17
Pt 18
Au 19
Hg 20
Tl 21
Pb 22
Bi 23
Po 24
At 25
Rn 26
EOF

cat Nval.txt SP.molden.input > SP.molden

### Calculate RESP charges
echo Running Multiwfn_noGUI...
Multiwfn_noGUI SP.molden -ispecial 1 > /dev/null << EOF
7
18
1
y
0
0
q
EOF

### Finalize
chgname=${1//$suffix/chg}
mv SP.chg $chgname
rm SP.* SP_* Nval.txt

echo Finished! The optimized atomic coordinates with RESP charges \(the last column\) have been exported to $chgname in current folder
echo Please properly cite Multiwfn_noGUI in your publication according to \"How to cite Multiwfn.pdf\" in Multiwfn package