// ENG1_QM_MOPAC.CPP

// Copyright (C) 2001 Tommi Hassinen.

// This package is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.

// This package is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.

// You should have received a copy of the GNU General Public License
// along with this package; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

/*################################################################################################*/

#include "libghemicalconfig2.h"
#include "eng1_qm_mopac.h"

#include "local_i18n.h"

#ifdef ENABLE_MOPAC7

#include <cstring>
#include <iomanip>
#include <sstream>
using namespace std;

#include <mopac7/libmopac7.h>

/*################################################################################################*/

eng1_qm_mopac * eng1_qm_mopac::mopac_lock = NULL;

eng1_qm_mopac::eng1_qm_mopac(setup * p1, i32u p2, i32u mode) : engine(p1, p2), eng1_qm(p1, p2)
{
	if (mopac_lock != NULL)
	{
		// the main MOPAC locking mechanism is now in model::CreateDefaultEngine().
		// if the locking fails in this stage, take it as a serious error and shutdown...
		
		assertion_failed(__FILE__, __LINE__, "mopac_lock != NULL");
	}
	else mopac_lock = this;
	
	if (GetSetup()->GetModel()->GetConstD_count() > 0)
	{
		GetSetup()->GetModel()->WarningMessage(CONSTRAINTS_NOT_SUPPORTED);
	}
	
	char fn_FOR005[256] = "FOR005";
	if (getenv("FOR005") != NULL) strcpy(fn_FOR005, getenv("FOR005"));
	
	cout << _("writing MOPAC-input file ") << fn_FOR005 << endl;
	
	ofstream file;
	file.open(fn_FOR005, ios::out);
	
	file << "XYZ NOLOG ";
	file << "SCFCRT=0.000001 ";	// slightly less strict than default, but seems to stabilize geomopt.
	
	switch (mode)
	{
		case MOPAC_MINDO3:
		file << "MINDO ";
		break;
		
		case MOPAC_AM1:
		file << "AM1 ";
		break;
		
		case MOPAC_PM3:
		file << "PM3 ";
		break;
		
		// if we don't write anything here, the default MNDO method will be used...
		// if we don't write anything here, the default MNDO method will be used...
		// if we don't write anything here, the default MNDO method will be used...
	};
	
// one should use UHF for transition state searching???
// but it will break the energy level diagrams...
// and unfortunately is quite slow...

//file << "UHF ";
//file << "TRIPLET ";
//file << "PULAY ";

	file << "CHARGE=" << GetSetup()->GetModel()->GetQMTotalCharge() << " ";
	file << "GEO-OK MMOK ";
	file << endl;
	
	file << "an automatically generated MOPAC input file." << endl << endl;
	
// add one dummy atom, to avoid program crash with molecules that have 3 first atoms linearly arranged.
// using this dummy atom won't have any other effect, except it cheats a test at GETGEO.F around line 357.
// removing this dummy atom and the test at GETGEO.F seems to produce exactly similar results...

// this problem can be studied with a molecule CH2=C=CH2 that can be drawn in either "right" way (carbons
// are the first three atoms) or "wrong" way... produces slightly higher energies and different MOs.

// HOW TO AVOID THIS PROBLEM: draw the molecule in 1-3-2 order instead of 1-2-3 order!!!
// HOW TO AVOID THIS PROBLEM: draw the molecule in 1-3-2 order instead of 1-2-3 order!!!
// HOW TO AVOID THIS PROBLEM: draw the molecule in 1-3-2 order instead of 1-2-3 order!!!

////////////////////////////////////////////////////////////////////////////////
//file << "XX     +1.00 1   +1.00 1   +1.00 1" << endl;	// this line can be commented out...
////////////////////////////////////////////////////////////////////////////////

	atom ** atmtab = GetSetup()->GetQMAtoms();
	for (i32s index = 0;index < GetSetup()->GetQMAtomCount();index++)
	{
		const fGL * defcrd = atmtab[index]->GetCRD(0);
		const fGL cf = 10.0;	// conversion factor for [nm] -> [Å]
		
		file << atmtab[index]->el.GetSymbol() << "\t";
		file << setprecision(6) << setw(12) << (defcrd[0] * cf) << " 1 ";
		file << setprecision(6) << setw(12) << (defcrd[1] * cf) << " 1 ";
		file << setprecision(6) << setw(12) << (defcrd[2] * cf) << " 1 ";
		file << endl;
	}
	
	file << endl;	// add an empty line to mark end of this input???
	file.close();
	
	// the MOPAC output is now directed to console; that's because at the beginning of LM7START()
	// at m7MAIN.c the unit 6 file is not opened. to send output to a logfile, open the unit 6 file.
	
	lm7start_();
	lm7_ini_full_xyz();
}

eng1_qm_mopac::~eng1_qm_mopac(void)
{
	if (mopac_lock != this) return;		// LOCK!!!
	
	lm7stop_();
	
	char fn_FOR005[256] = "FOR005";
	if (getenv("FOR005") != NULL) strcpy(fn_FOR005, getenv("FOR005"));
	
	char fn_SHUTDOWN[256] = "SHUTDOWN";
	if (getenv("SHUTDOWN") != NULL) strcpy(fn_SHUTDOWN, getenv("SHUTDOWN"));
	
	cout << _("removing intermediate files... ");
	
	static ostringstream cs_FOR005;
	cs_FOR005 << "rm " << fn_FOR005 << ends;
	system(cs_FOR005.str().c_str());
	
	static ostringstream cs_SHUTDOWN;
	cs_SHUTDOWN << "rm " << fn_SHUTDOWN << ends;
	system(cs_SHUTDOWN.str().c_str());
	
	cout << "OK!" << endl;
	mopac_lock = NULL;
}

void eng1_qm_mopac::FixMeLaterTSS(void)
{
// this is a temporary fix for some VERY bad problems in TSS...
// this is a temporary fix for some VERY bad problems in TSS...
// this is a temporary fix for some VERY bad problems in TSS...
	
	lm7stop_();
	
	char fn_FOR005[256] = "FOR005";
	if (getenv("FOR005") != NULL) strcpy(fn_FOR005, getenv("FOR005"));
	char fn_SHUTDOWN[256] = "SHUTDOWN";
	if (getenv("SHUTDOWN") != NULL) strcpy(fn_SHUTDOWN, getenv("SHUTDOWN"));
	
	static ostringstream cs_FOR005;
	cs_FOR005 << "rm " << fn_FOR005 << ends;
	system(cs_FOR005.str().c_str());
	static ostringstream cs_SHUTDOWN;
	cs_SHUTDOWN << "rm " << fn_SHUTDOWN << ends;
	system(cs_SHUTDOWN.str().c_str());
	
	mopac_lock = NULL;
}

i32s eng1_qm_mopac::GetElectronCount(void)
{
	return lm7_get_electron_count();
}

i32s eng1_qm_mopac::GetOrbitalCount(void)
{
	return lm7_get_orbital_count();
}

f64 eng1_qm_mopac::GetOrbitalEnergy(i32s orb_i)
{
	return lm7_get_orbital_energy(orb_i);
}

void eng1_qm_mopac::Compute(i32u p1, bool)
{
	if (mopac_lock != this) return;		// LOCK!!!
	
	double e; double hf; int i;
	
	for (i = 0;i < lm7_get_atom_count();i++)
	{
		lm7_set_atom_crd(i, & crd[l2g_qm[i] * 3]);
	}
	
	if (p1 == 0)		// energy was requested...
	{
		lm7_call_compfg(& e, & hf, 0);
	}
	else if (p1 == 1)	// energy and gradient was requested...
	{
		lm7_call_compfg(& e, & hf, 1);
		
		for (i = 0;i < lm7_get_atom_count();i++)
		{
			lm7_get_atom_grad(i, & d1[l2g_qm[i] * 3]);
		}
	}
	else	// can't calculate higher derivatives just yet...
	{
		assertion_failed(__FILE__, __LINE__, "not_implemented");
	}
	
//cout << "energy = " << e << " kJ/mol" << endl;
//cout << "heat of formation = " << hf << " kcal/mol" << endl;

	energy = e;
	
	// the rest is for transition state search only!!!
	// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	
	if (tss_ref_str == NULL) return;
	
	// this is currently implemented for MOPAC only; there is no problems in
	// making it for MPQC as well, but since semi-empirical methods are quicker
	// this seems to be enough. ONE CAN USE AB INITIO IN STATIONARY STATE SEARCH!
	
	tss_delta_ene = 0.0;
	for (i = 0;i < lm7_get_atom_count();i++)
	{
		f64 t1a[3]; f64 t1b = 0.0;
		for (i32s n5 = 0;n5 < 3;n5++)
		{
			f64 t9a = crd[l2g_qm[i] * 3 + n5];
			f64 t9b = tss_ref_str[i * 3 + n5];
			
			t1a[n5] = t9a - t9b;
			t1b += t1a[n5] * t1a[n5];
		}
		
		f64 t1c = sqrt(t1b);
		for (i32s n5 = 0;n5 < 3;n5++)
		{
			f64 t9a = t1a[n5] / t1c;
			
			t1a[n5] = +t9a;
		}
		
		// f = a(x)^2
		// df/dx = 2a(x)
		
		f64 t2a = tss_force_const * t1b;
		energy += t2a;
		
		// df/fa = (x)^2
		
		tss_delta_ene += t1b;
		
		if (p1 > 0)
		{
			f64 t2b = 2.0 * tss_force_const * t1c;
			
			for (i32s n5 = 0;n5 < 3;n5++)
			{
				f64 t2c = t1a[n5] * t2b;
				
			//	d1[i * 3 + n5] += t2c;
				d1[l2g_qm[i] * 3 + n5] += t2c;
			}
		}
	}
}

void eng1_qm_mopac::SetupPlotting(void)
{
	lm7iniplt_();
}

fGL eng1_qm_mopac::GetESP(fGL * pp, fGL * dd)
{
	if (mopac_lock != this) return 0.0;	// LOCK!!!
	
	const double cf1 = 10.0;	// conversion factor for [nm] -> [Å]
	double tmpcrd[3] = { pp[0] * cf1, pp[1] * cf1, pp[2] * cf1 };
	
	lm7_set_num_potesp(1);
	lm7_set_crd_potesp(0, tmpcrd);
	
	getesp_();
	
	const double cf2 = 2625.5;	// conversion factor for [Hartree] -> [kJ/mol]
	fGL value = lm7_get_val_potesp(0) * cf2;
	
	if (dd != NULL)		// numerical gradient...
	{
		fGL old;
		const fGL delta = 0.0001;
		
		old = pp[0]; pp[0] += delta;
		dd[0] = (GetESP(pp, NULL) - value) / delta;
		pp[0] = old;
		
		old = pp[1]; pp[1] += delta;
		dd[1] = (GetESP(pp, NULL) - value) / delta;
		pp[1] = old;
		
		old = pp[2]; pp[2] += delta;
		dd[2] = (GetESP(pp, NULL) - value) / delta;
		pp[2] = old;
	}
	
	return value;
}

fGL eng1_qm_mopac::GetElDens(fGL * pp, fGL * dd)
{
	if (mopac_lock != this) return 0.0;	// LOCK!!!
	
	const double cf = 18.897162;	// conversion factor for [nm] -> [bohr]
	double tmpcrd[3] = { pp[0] * cf, pp[1] * cf, pp[2] * cf };
	
	lm7_set_num_potesp(1);
	lm7_set_crd_potesp(0, tmpcrd);
	
	geteldens_();
	
	fGL value = lm7_get_val_potesp(0);
	
	if (dd != NULL)		// numerical gradient...
	{
		fGL old;
		const fGL delta = 0.0001;
		
		old = pp[0]; pp[0] += delta;
		dd[0] = (GetElDens(pp, NULL) - value) / delta;
		pp[0] = old;
		
		old = pp[1]; pp[1] += delta;
		dd[1] = (GetElDens(pp, NULL) - value) / delta;
		pp[1] = old;
		
		old = pp[2]; pp[2] += delta;
		dd[2] = (GetElDens(pp, NULL) - value) / delta;
		pp[2] = old;
	}
	
	return value;
}

fGL eng1_qm_mopac::GetOrbital(fGL * pp, fGL * dd)
{
	if (mopac_lock != this) return 0.0;	// LOCK!!!
	
	const double cf = 18.897162;	// conversion factor for [nm] -> [bohr]
	double tmpcrd[3] = { pp[0] * cf, pp[1] * cf, pp[2] * cf };
	
	lm7_set_num_potesp(1);
	lm7_set_crd_potesp(0, tmpcrd);
	
	lm7_set_plots_orbital_index(GetSetup()->GetModel()->qm_current_orbital + 1);
	
	getorb_();
	
	fGL value = lm7_get_val_potesp(0);
	
	if (dd != NULL)		// numerical gradient...
	{
		fGL old;
		const fGL delta = 0.0001;
		
		old = pp[0]; pp[0] += delta;
		dd[0] = (GetOrbital(pp, NULL) - value) / delta;
		pp[0] = old;
		
		old = pp[1]; pp[1] += delta;
		dd[1] = (GetOrbital(pp, NULL) - value) / delta;
		pp[1] = old;
		
		old = pp[2]; pp[2] += delta;
		dd[2] = (GetOrbital(pp, NULL) - value) / delta;
		pp[2] = old;
	}
	
	return value;
}

fGL eng1_qm_mopac::GetOrbDens(fGL * pp, fGL * dd)
{
	if (mopac_lock != this) return 0.0;	// LOCK!!!
	
	// this orbital density plot is closely related to the orbital plot;
	// we just square the orbital function and multiply it with 2 (assuming 2 electrons per orbital, true for RHF).
	
	fGL orb = GetOrbital(pp, dd);
	fGL value = 2.0 * orb * orb;
	
	if (dd != NULL)		// numerical gradient...
	{
		fGL old;
		const fGL delta = 0.0001;
		
		old = pp[0]; pp[0] += delta;
		dd[0] = (GetOrbDens(pp, NULL) - value) / delta;
		pp[0] = old;
		
		old = pp[1]; pp[1] += delta;
		dd[1] = (GetOrbDens(pp, NULL) - value) / delta;
		pp[1] = old;
		
		old = pp[2]; pp[2] += delta;
		dd[2] = (GetOrbDens(pp, NULL) - value) / delta;
		pp[2] = old;
	}
	
	return value;
}

/*################################################################################################*/

#endif	// ENABLE_MOPAC7

// eof