## Getting ready
## =============
cd ..
cp -a tutorial.orig tutorial
cd tutorial
bash
export CONPATH=`pwd`/scripts
PATH=$PATH:$CONPATH

## TDAZ
## ====
# get toppar/tdaz_init.str from paramchem.org
cd models/tdaz
# Generate Gaussian input
charmm str:../../toppar/tdaz_init.str -i generate.inp -o generate.out

# Run file created in gauss subdirectory with Gaussian
# Use resulting optimized geometry to start water interaction

# Constrained charge optimization
charmm str:../../toppar/tdaz_init.str -i water_constr.inp -o water_constr_init.out
cp -p tdaz_water_constr.ene tdaz_water_constr_init.ene
cat tdaz_water_constr_init.ene
charmm str:../../toppar/tdaz_water1.str -i water_constr.inp -o water_constr_min.out
cat tdaz_water_constr.ene
cp -p tdaz_water_constr.ene tdaz_water_constr_min.ene

# Quick
charmm str:../../toppar/tdaz_water1.str -i quick_molvib.inp -o quick_molvib_init.out
less +/quick quick_molvib_init.out
angsum ../../toppar/par_all36_cgenff.prm ../../toppar/tdaz_water1.str 0 TDAZ ../../toppar/angsum1.rtf
# RESI TDAZ ring 5 C1 N2 N3 C4 S5 : sum = 545 (540)
# Parameters involved:
# ../../toppar/tdaz_water1.str line 37: CG2R53 NG2R50 NG2R50 ang=106.80
# ../../toppar/tdaz_water1.str line 37: CG2R53 NG2R50 NG2R50 ang=106.80
# ../../toppar/par_all36_cgenff.prm line 1059: NG2R50 CG2R53 SG2R50 ang=117.20
# ../../toppar/tdaz_water1.str line 38: CG2R53 SG2R50 CG2R53 ang=97.00
# ../../toppar/par_all36_cgenff.prm line 1059: NG2R50 CG2R53 SG2R50 ang=117.20
angsum ../../toppar/par_all36_cgenff.prm ../../toppar/tdaz_quick1.str 0 TDAZ ../../toppar/angsum1.rtf
charmm str:../../toppar/tdaz_quick1.str -i quick_molvib.inp -o quick_molvib_quick1.out
less +/quick quick_molvib_quick1.out

# Molvib
charmm -i tdaz_molvib_qm.inp -o tdaz_molvib_qm.out
getmolvib tdaz_molvib_qm.out
getmolvib quick_molvib_quick1.out
charmm str:../../toppar/tdaz_quick1molvib.str -i quick_molvib.inp -o quick_molvib_quick1molvib.out
getmolvib quick_molvib_quick1molvib.out

# Revisit quick
less +/quick quick_molvib_quick1molvib.out
charmm str:../../toppar/tdaz_quick2molvib.str -i quick_molvib.inp -o quick_molvib_opt.out
less +/quick quick_molvib_opt.out
getmolvib quick_molvib_opt.out
# self-consistent

# Relaxed charge optimization
charmm str:../../toppar/tdaz_quick1molvib.str -i water.inp -o water.out
cat tdaz_water_constr_min.ene
cat tdaz_water.ene
# self-consistent
cd ../..

## SULFAZ
## ======
# manually put together toppar/acazam_models_init.str

cd models/sulfaz
mkdir -p fit_ang fit_dih

# Prepare gaussian outputs
charmm str:../../toppar/acazam_models_init.str -i generate_dummy.inp -o generate_dummy.out
cd scanang
awk '/Scan/{print $3; exit}' sulfaz_scanang.log # shortcut to get the scan string for the following line
pescrd sulfaz_scanang.log ../sulfaz_dummy.crd EUMP2 'A(4,5,7)'
echo >> s_file.txt
cp -p e_file.txt ../fit_ang/scan_sulfaz.qme
cd ../scandih
awk '/Scan/{print $3; exit}' sulfaz_scan.log
pescrd sulfaz_scan.log ../sulfaz_dummy.crd EUMP2 'D(6,5,7,10)'
echo >> s_file.txt
cp -p e_file.txt ../fit_dih/scan_sulfaz.qme
cd ..

# Fit angle
cd fit_ang
./initialize
charmm str:../../../toppar/acazam_models_init.str quick:1 -i angscan.inp -o angscan_init.out
grep exceeded angscan_init.out
## The following is normal, but anything else not:
# CONJUG> Minimization exiting with number of steps limit (  200) exceeded.
grep 'A(4,5,7)' ../gauss/sulfaz_opt_freq_mp2b.log
# 121.8386
cat scan_sulfaz.quick
# 121.093
angsum ../../../toppar/par_all36_cgenff.prm ../../../toppar/acazam_models_init.str 0 SULFAZ ../../../toppar/angsum2.rtf
# RESI SULFAZ center C4 bound to N3 S5 S9 : sum = 335.2 (360)
# Parameters involved:
# ../../../toppar/par_all36_cgenff.prm line 1059: NG2R50 CG2R53 SG2R50 ang=117.20
# ../../../toppar/acazam_models_init.str line 145: SG2R50 CG2R53 SG3O2  ang=109.00
# ../../../toppar/acazam_models_init.str line 143: NG2R50 CG2R53 SG3O2  ang=109.00
# RESI SULFAZ center N10 bound to S9 H101 H102 : sum = 340 (360)
# Parameters involved:
# ../../../toppar/par_all36_cgenff.prm line 1924: SG3O2  NG321  HGP1   ang=115.00
# ../../../toppar/par_all36_cgenff.prm line 1925: HGP1   NG321  HGP1   ang=110.00
# ../../../toppar/par_all36_cgenff.prm line 1924: SG3O2  NG321  HGP1   ang=115.00
angsum ../../../toppar/par_all36_cgenff.prm ../../../toppar/acazam_sulfaz_quick1.str 0 SULFAZ ../../../toppar/angsum2.rtf
# RESI SULFAZ center N10 bound to S9 H101 H102 : sum = 340 (360)
# Parameters involved:
# ../../../toppar/par_all36_cgenff.prm line 1924: SG3O2  NG321  HGP1   ang=115.00
# ../../../toppar/par_all36_cgenff.prm line 1925: HGP1   NG321  HGP1   ang=110.00
# ../../../toppar/par_all36_cgenff.prm line 1924: SG3O2  NG321  HGP1   ang=115.00
./initialize
charmm str:../../../toppar/acazam_sulfaz_quick1.str -i angscan.inp -o angscan_quick1.out
cat scan_sulfaz.quick
# 121.767
scanplot -i init scan_sulfaz.qme scan_sulfaz.mme
./initialize
charmm str:../../../toppar/acazam_sulfaz_quick1fc.str -i angscan.inp -o angscan_quick1fc.out
scanplot -i quick1fc scan_sulfaz.qme scan_sulfaz.mme
cat scan_sulfaz.quick
# 122.209
./initialize
charmm str:../../../toppar/acazam_sulfaz_quick2fc.str -i angscan.inp -o angscan_quick2fc.out
cat scan_sulfaz.quick
# 121.858
scanplot -i quick2fc scan_sulfaz.qme scan_sulfaz.mme
# self-consistent
cd ..

# Fit dihedral
cd fit_dih
./initialize
charmm str:../../../toppar/acazam_sulfaz_quick2fc.str -i allscan.inp -o allscan_init.out
grep 'NRAPH.\+exceeded' allscan_init.out
scanplot -i init scan_sulfaz.qme scan_sulfaz.mme
./initialize
charmm str:../../../toppar/acazam_sulfaz_3fold.str -i allscan.inp -o allscan_3fold.out
scanplot -i 3fold scan_sulfaz.qme scan_sulfaz.mme
./initialize
charmm str:../../../toppar/acazam_sulfaz_1fold.str -i allscan.inp -o allscan_1fold.out
scanplot -i 1fold scan_sulfaz.qme scan_sulfaz.mme
# Rebalanced 3-fold and 1-fold and introduced 6-fold
./initialize
charmm str:../../../toppar/acazam_sulfaz_6fold.str -i allscan.inp -o allscan_6fold.out
scanplot -i 6fold scan_sulfaz.qme scan_sulfaz.mme
# Used O. Guvench' MCSA fitting to introduce 4-fold that shifts location of minimum
./initialize
charmm str:../../../toppar/acazam_sulfaz_4fold.str -i allscan.inp -o allscan_4fold.out
scanplot -i 4fold scan_sulfaz.qme scan_sulfaz.mme
cd ../../..

## TDAZAM
## ======
cd models/tdazam
mkdir -p fit_ang fit_dih fit_impr

# Prepare gaussian outputs
charmm str:../../toppar/acazam_sulfaz_4fold.str -i generate_dummy.inp -o generate_dummy.out
cd scanang
awk '/Scan/{print $3; exit}' tdazam_scanang2.log
pescrd tdazam_scanang2.log ../tdazam_dummy.crd EUMP2 'A(2,1,7)'
echo >> s_file.txt
cp -p e_file.txt ../fit_ang/scan_tdazam.qme
cd ../scanimp
awk '/Scan/{print $3; exit}' tdazam_scanimp.log
pescrd tdazam_scanimp.log ../tdazam_dummy.crd EUMP2 'D(1,2,6,7)'
echo >> s_file.txt
cp -p e_file.txt ../fit_imp/scan_tdazam.qme
cd ../scandih
awk '/Scan/{print $3; exit}' tdazam_scan.log
pescrd tdazam_scan.log ../tdazam_dummy.crd EUMP2 'D(2,1,7,9)'
echo >> s_file.txt
cp -p e_file.txt ../fit_dih/scan_tdazam.qme
cd ..

# Fit angle
cd fit_ang
./initialize
charmm str:../../../toppar/acazam_sulfaz_4fold.str quick:1 -i angscan.inp -o angscan_init.out
grep exceeded angscan_init.out
## The following is normal, but anything else not:
# CONJUG> Minimization exiting with number of steps limit (  200) exceeded.
grep 'A(2,1,7)' ../scanimp/tdazam_scanimp.log
# 119.7428
## Sorry, accidentally deleted a file. If I hadn't, the above would have been:
## grep 'A(2,1,7)' ../gauss/tdazam_opt_freq_mp2.log
## 119.7494
cat scan_tdazam.quick
# 117.955
angsum ../../../toppar/par_all36_cgenff.prm ../../../toppar/acazam_sulfaz_4fold.str 0 TDAZAM ../../../toppar/angsum2.rtf
# RESI TDAZAM center C1 bound to S5 N2 N6 : sum = 363 (360)
# Parameters involved:
# ../../../toppar/par_all36_cgenff.prm line 1059: NG2R50 CG2R53 SG2R50 ang=117.20
# ../../../toppar/acazam_sulfaz_4fold.str line 145: NG2R50 CG2R53 NG2S1  ang=127.80
# ../../../toppar/acazam_sulfaz_4fold.str line 146: NG2S1  CG2R53 SG2R50 ang=118.00
angsum ../../../toppar/par_all36_cgenff.prm ../../../toppar/acazam_tdazam_quick1.str 0 TDAZAM ../../../toppar/angsum2.rtf
./initialize
charmm str:../../../toppar/acazam_tdazam_quick1.str -i angscan.inp -o angscan_quick1.out
cat scan_tdazam.quick
# 119.738
scanplot -i init scan_tdazam.qme scan_tdazam.mme
./initialize
charmm str:../../../toppar/acazam_tdazam_quick1fc.str -i angscan.inp -o angscan_quick1fc.out
scanplot -i quick1fc scan_tdazam.qme scan_tdazam.mme
cat scan_tdazam.quick
# 122.355
./initialize
charmm str:../../../toppar/acazam_tdazam_quick2fc.str -i angscan.inp -o angscan_quick2fc.out
cat scan_tdazam.quick
# 119.776
scanplot -i quick2fc scan_tdazam.qme scan_tdazam.mme
# self-consistent
cd ..

# Fit improper
cd fit_imp
./initialize
charmm str:../../../toppar/acazam_tdazam_quick2fc.str -i allscan.inp -o allscan_init.out
grep 'NRAPH.\+exceeded' allscan_init.out
scanplot -i init scan_tdazam.qme scan_tdazam.mme
./initialize
charmm str:../../../toppar/acazam_tdazam_imp.str -i allscan.inp -o allscan_opt.out
scanplot -i opt scan_tdazam.qme scan_tdazam.mme
cd ..

# Fit dihedral
cd fit_dih
./initialize
charmm str:../../../toppar/acazam_tdazam_imp.str -i allscan.inp -o allscan_init.out
grep 'NRAPH.\+exceeded' allscan_init.out
scanplot -i init scan_tdazam.qme scan_tdazam.mme
./initialize
charmm str:../../../toppar/acazam_tdazam_1fold.str -i allscan.inp -o allscan_1fold.out
scanplot -i 1fold scan_tdazam.qme scan_tdazam.mme
./initialize
charmm str:../../../toppar/acazam_tdazam_2fold.str -i allscan.inp -o allscan_2fold.out
scanplot -i 2fold scan_tdazam.qme scan_tdazam.mme
./initialize
charmm str:../../../toppar/acazam_tdazam_3fold.str -i allscan.inp -o allscan_3fold.out
scanplot -i 3fold scan_tdazam.qme scan_tdazam.mme
./initialize
charmm str:../../../toppar/acazam_tdazam_4fold.str -i allscan.inp -o allscan_4fold.out
scanplot -i 4fold scan_tdazam.qme scan_tdazam.mme
cd ../../..
