OSDN -- Open-Source-Software-Entwicklung und Downloads
Übersetzungsstatus von Deutsch translation status for de
  • R/O
  • HTTP
  • SSH
  • HTTPS
Fork

Commit

Tags
Keine Tags

Frequently used words (click to add to your profile)

javac++androidlinuxc#windowsobjective-ccocoa誰得qtpythonphprubygameguibathyscaphec計画中(planning stage)翻訳omegatframeworktwitterdomtestvb.netdirectxゲームエンジンbtronarduinopreviewer

Commit MetaInfo

Revision2a09f969ccaf8f907b0c7f1491f40b741874a61f (tree)
Zeit2013-11-29 16:51:34
AutorMikiya Fujii <mikiya.fujii@gmai...>
CommiterMikiya Fujii

Log Message

EnvironmentalPointCharge-class(EPC) is added. In molecule-class distance matrix for EPC is added. InputParser-class is modified for EPCs. Some methods for two electrons integral is refactored such as rename. README is updated. #32469

git-svn-id: https://svn.sourceforge.jp/svnroot/molds/trunk@1569 1136aad2-a195-0410-b898-f5ea1d11b9d8

Ändern Zusammenfassung

Diff

--- a/doc/README.txt
+++ b/doc/README.txt
@@ -338,6 +338,43 @@ HOW TO WRITE INPUT:
338338 file_prefix MOPlot_
339339 MOPLOT_END
340340
341+ <Environmental Point Charge(EPC) method>
342+ Environmental point charge method is a simplified method of the QM/MM
343+ because the environmental point changes are treated as atoms in the MM region.
344+ The differences between the QM/MM and EPC are summarized below:
345+ - Electrostatic interaction between QM and MM region:
346+ QM/MM: Electrostatic interaction is mutually added to QM and MM atoms.
347+ EPC : Electrostatic field caused by the EPCs affects the QM region
348+ although the each EPC is not affected by electrostatic field
349+ caused by the QM atoms and other EPCs.
350+ Namely, each EPC is fixed at point of space.
351+ - Van der Waals interaction between QM and MM region:
352+ QM/MM: Included.
353+ EPC : Not included.
354+ This EPC method can be used with MNDO-series (MNDO, AM1, AM1-D, PM3, PM3-D, and PDDG/PM3) only.
355+ To use this environmental point charges method, write EPC-directive.
356+
357+ E.g.
358+ EPC
359+ (options)
360+ EPC_END
361+
362+ -options
363+ "the cartesian coordinates and charge" is only prepared.
364+ Namely, each line should containe 4 doubles.
365+ The first three doubles are the cartesian coordinates of
366+ each environmental point charge in angstrom unit.
367+ The last double is the charge in atomic unit,
368+ e.g. -1 and 1 mean charge of an electron and a proton, respectively.
369+ Multiple setting of the environmental point charge is approvable, of course.
370+
371+ E.g.
372+ EPC
373+ 0.0 0.0 0.0 -1
374+ 2.2 1.5 3.0 -1.5
375+ 0.0 2.0 5.0 0.3
376+ EPC_END
377+
341378 <Frequencies (Normal modes analysis)>
342379 write frequencies-directive. Note taht not only the frequencies but also the normal modes are calculated.
343380
--- a/src/Makefile
+++ b/src/Makefile
@@ -37,9 +37,9 @@ EXENAME = molds
3737 DEPFILE = obj/objfile.dep
3838 LDFLAGS =
3939
40-ALL_CPP_FILES = base/Enums.cpp base/PrintController.cpp base/MolDSException.cpp base/MallocerFreer.cpp mpi/MpiProcess.cpp mpi/AsyncCommunicator.cpp wrappers/Blas.cpp wrappers/Lapack.cpp base/Utilities.cpp base/MathUtilities.cpp base/EularAngle.cpp base/Parameters.cpp base/RealSphericalHarmonicsIndex.cpp base/atoms/Atom.cpp base/atoms/Hatom.cpp base/atoms/Liatom.cpp base/atoms/Catom.cpp base/atoms/Natom.cpp base/atoms/Oatom.cpp base/atoms/Satom.cpp base/factories/AtomFactory.cpp base/Molecule.cpp base/InputParser.cpp base/GTOExpansionSTO.cpp base/loggers/MOLogger.cpp base/loggers/DensityLogger.cpp base/loggers/HoleDensityLogger.cpp base/loggers/ParticleDensityLogger.cpp cndo/Cndo2.cpp indo/Indo.cpp zindo/ZindoS.cpp mndo/Mndo.cpp am1/Am1.cpp am1/Am1D.cpp pm3/Pm3.cpp pm3/Pm3D.cpp pm3/Pm3Pddg.cpp base/factories/ElectronicStructureFactory.cpp md/MD.cpp mc/MC.cpp rpmd/RPMD.cpp nasco/NASCO.cpp optimization/Optimizer.cpp optimization/ConjugateGradient.cpp optimization/SteepestDescent.cpp optimization/BFGS.cpp optimization/GEDIIS.cpp base/factories/OptimizerFactory.cpp base/MolDS.cpp Main.cpp
41-ALL_HEAD_FILES = config.h base/Enums.h base/Uncopyable.h base/PrintController.h base/MolDSException.h base/containers/ThreadSafeQueue.h base/MallocerFreer.h mpi/MpiInt.h mpi/MpiProcess.h mpi/AsyncCommunicator.h wrappers/Blas.h wrappers/Lapack.h base/Utilities.h base/MathUtilities.h base/EularAngle.h base/Parameters.h base/RealSphericalHarmonicsIndex.h base/atoms/Atom.h base/atoms/Hatom.h base/atoms/Liatom.h base/atoms/Catom.h base/atoms/Natom.h base/atoms/Oatom.h base/atoms/Satom.h base/factories/AtomFactory.h base/Molecule.h base/InputParser.h base/GTOExpansionSTO.h base/loggers/MOLogger.h base/loggers/DensityLogger.h base/loggers/HoleDensityLogger.h base/loggers/ParticleDensityLogger.h base/ElectronicStructure.h cndo/Cndo2.h cndo/ReducedOverlapAOsParameters.h indo/Indo.h zindo/ZindoS.h mndo/Mndo.h am1/Am1.h am1/Am1D.h pm3/Pm3.h pm3/Pm3D.h pm3/Pm3Pddg.h base/factories/ElectronicStructureFactory.h md/MD.h mc/MC.h rpmd/RPMD.h nasco/NASCO.h optimization/Optimizer.h optimization/ConjugateGradient.h optimization/SteepestDescent.h optimization/BFGS.h optimization/GEDIIS.h base/factories/OptimizerFactory.h base/MolDS.h
42-ALL_OBJ_FILES = obj/Enums.o obj/PrintController.o obj/MolDSException.o obj/MallocerFreer.o obj/MpiProcess.o obj/AsyncCommunicator.o obj/Blas.o obj/Lapack.o obj/Utilities.o obj/MathUtilities.o obj/EularAngle.o obj/Parameters.o obj/RealSphericalHarmonicsIndex.o obj/Atom.o obj/Hatom.o obj/Liatom.o obj/Catom.o obj/Natom.o obj/Oatom.o obj/Satom.o obj/AtomFactory.o obj/Molecule.o obj/InputParser.o obj/GTOExpansionSTO.o obj/MOLogger.o obj/DensityLogger.o obj/HoleDensityLogger.o obj/ParticleDensityLogger.o obj/Cndo2.o obj/Indo.o obj/ZindoS.o obj/Mndo.o obj/Am1.o obj/Am1D.o obj/Pm3.o obj/Pm3D.o obj/Pm3Pddg.o obj/ElectronicStructureFactory.o obj/MD.o obj/MC.o obj/RPMD.o obj/NASCO.o obj/Optimizer.o obj/ConjugateGradient.o obj/SteepestDescent.o obj/BFGS.o obj/GEDIIS.o obj/OptimizerFactory.o obj/MolDS.o obj/Main.o
40+ALL_CPP_FILES = base/Enums.cpp base/PrintController.cpp base/MolDSException.cpp base/MallocerFreer.cpp mpi/MpiProcess.cpp mpi/AsyncCommunicator.cpp wrappers/Blas.cpp wrappers/Lapack.cpp base/Utilities.cpp base/MathUtilities.cpp base/EularAngle.cpp base/Parameters.cpp base/RealSphericalHarmonicsIndex.cpp base/atoms/Atom.cpp base/atoms/Hatom.cpp base/atoms/Liatom.cpp base/atoms/Catom.cpp base/atoms/Natom.cpp base/atoms/Oatom.cpp base/atoms/Satom.cpp base/atoms/mm/EnvironmentalPointCharge.cpp base/factories/AtomFactory.cpp base/Molecule.cpp base/InputParser.cpp base/GTOExpansionSTO.cpp base/loggers/MOLogger.cpp base/loggers/DensityLogger.cpp base/loggers/HoleDensityLogger.cpp base/loggers/ParticleDensityLogger.cpp cndo/Cndo2.cpp indo/Indo.cpp zindo/ZindoS.cpp mndo/Mndo.cpp am1/Am1.cpp am1/Am1D.cpp pm3/Pm3.cpp pm3/Pm3D.cpp pm3/Pm3Pddg.cpp base/factories/ElectronicStructureFactory.cpp md/MD.cpp mc/MC.cpp rpmd/RPMD.cpp nasco/NASCO.cpp optimization/Optimizer.cpp optimization/ConjugateGradient.cpp optimization/SteepestDescent.cpp optimization/BFGS.cpp optimization/GEDIIS.cpp base/factories/OptimizerFactory.cpp base/MolDS.cpp Main.cpp
41+ALL_HEAD_FILES = config.h base/Enums.h base/Uncopyable.h base/PrintController.h base/MolDSException.h base/containers/ThreadSafeQueue.h base/MallocerFreer.h mpi/MpiInt.h mpi/MpiProcess.h mpi/AsyncCommunicator.h wrappers/Blas.h wrappers/Lapack.h base/Utilities.h base/MathUtilities.h base/EularAngle.h base/Parameters.h base/RealSphericalHarmonicsIndex.h base/atoms/Atom.h base/atoms/Hatom.h base/atoms/Liatom.h base/atoms/Catom.h base/atoms/Natom.h base/atoms/Oatom.h base/atoms/Satom.h base/atoms/mm/EnvironmentalPointCharge.h base/factories/AtomFactory.h base/Molecule.h base/InputParser.h base/GTOExpansionSTO.h base/loggers/MOLogger.h base/loggers/DensityLogger.h base/loggers/HoleDensityLogger.h base/loggers/ParticleDensityLogger.h base/ElectronicStructure.h cndo/Cndo2.h cndo/ReducedOverlapAOsParameters.h indo/Indo.h zindo/ZindoS.h mndo/Mndo.h am1/Am1.h am1/Am1D.h pm3/Pm3.h pm3/Pm3D.h pm3/Pm3Pddg.h base/factories/ElectronicStructureFactory.h md/MD.h mc/MC.h rpmd/RPMD.h nasco/NASCO.h optimization/Optimizer.h optimization/ConjugateGradient.h optimization/SteepestDescent.h optimization/BFGS.h optimization/GEDIIS.h base/factories/OptimizerFactory.h base/MolDS.h
42+ALL_OBJ_FILES = obj/Enums.o obj/PrintController.o obj/MolDSException.o obj/MallocerFreer.o obj/MpiProcess.o obj/AsyncCommunicator.o obj/Blas.o obj/Lapack.o obj/Utilities.o obj/MathUtilities.o obj/EularAngle.o obj/Parameters.o obj/RealSphericalHarmonicsIndex.o obj/Atom.o obj/Hatom.o obj/Liatom.o obj/Catom.o obj/Natom.o obj/Oatom.o obj/Satom.o obj/EnvironmentalPointCharge.o obj/AtomFactory.o obj/Molecule.o obj/InputParser.o obj/GTOExpansionSTO.o obj/MOLogger.o obj/DensityLogger.o obj/HoleDensityLogger.o obj/ParticleDensityLogger.o obj/Cndo2.o obj/Indo.o obj/ZindoS.o obj/Mndo.o obj/Am1.o obj/Am1D.o obj/Pm3.o obj/Pm3D.o obj/Pm3Pddg.o obj/ElectronicStructureFactory.o obj/MD.o obj/MC.o obj/RPMD.o obj/NASCO.o obj/Optimizer.o obj/ConjugateGradient.o obj/SteepestDescent.o obj/BFGS.o obj/GEDIIS.o obj/OptimizerFactory.o obj/MolDS.o obj/Main.o
4343
4444 $(EXENAME): $(ALL_OBJ_FILES)
4545 $(CC) -o $@ -Wl,-rpath=$(BOOST_LIB_DIR) -L$(BOOST_LIB_DIR) $(LDFLAGS) $(ALL_OBJ_FILES) $(LIBS)
--- a/src/am1/Am1.cpp
+++ b/src/am1/Am1.cpp
@@ -92,20 +92,20 @@ void Am1::SetMessages(){
9292 = "Error in am1::Am1::GetNddoRepulsionIntegral1stDerivative: Bad orbital is set.\n";
9393 this->errorMessageGetNddoRepulsionIntegral2ndDerivative
9494 = "Error in am1::Am1::GetNddoRepulsionIntegral2ndDerivative: Bad orbital is set.\n";
95- this->errorMessageCalcTwoElecTwoCoreNullMatrix
96- = "Error in am1::Am1::CalcTwoElecTwoCore: The two elec two core matrix is NULL.\n";
97- this->errorMessageCalcDiatomicTwoElecTwoCoreSameAtoms
98- = "Error in am1::Am1::CalcDiatomicTwoElecTwoCore: Atom A and B is same.\n";
99- this->errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesSameAtoms
100- = "Error in am1::Am1::CalcDiatomicTwoElecTwoCore1stDerivatives: Atom A and B is same.\n";
101- this->errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesSameAtoms
102- = "Error in am1::Am1::CalcDiatomicTwoElecTwoCore2ndDerivatives: Atom A and B is same.\n";
103- this->errorMessageCalcDiatomicTwoElecTwoCoreNullMatrix
104- = "Error in am1::Am1::CalcDiatomicTwoElecTwoCore: The two elec two core diatomic matrix is NULL.\n";
105- this->errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesNullMatrix
106- = "Error in am1::Am1::CalcDiatomicTwoElecTwoCore1stDerivatives: The two elec two core diatomic matrix is NULL.\n";
107- this->errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesNullMatrix
108- = "Error in am1::Am1::CalcDiatomicTwoElecTwoCore2ndDerivatives: The two elec two core diatomic matrix is NULL.\n";
95+ this->errorMessageCalcTwoElecsTwoCoresNullMatrix
96+ = "Error in am1::Am1::CalcTwoElecsTwoCores: The two elec two core matrix is NULL.\n";
97+ this->errorMessageCalcDiatomicTwoElecsTwoCoresSameAtoms
98+ = "Error in am1::Am1::CalcDiatomicTwoElecsTwoCores: Atom A and B is same.\n";
99+ this->errorMessageCalcDiatomicTwoElecsTwoCores1stDerivativesSameAtoms
100+ = "Error in am1::Am1::CalcDiatomicTwoElecsTwoCores1stDerivatives: Atom A and B is same.\n";
101+ this->errorMessageCalcDiatomicTwoElecsTwoCores2ndDerivativesSameAtoms
102+ = "Error in am1::Am1::CalcDiatomicTwoElecsTwoCores2ndDerivatives: Atom A and B is same.\n";
103+ this->errorMessageCalcDiatomicTwoElecsTwoCoresNullMatrix
104+ = "Error in am1::Am1::CalcDiatomicTwoElecsTwoCores: The two elec two core diatomic matrix is NULL.\n";
105+ this->errorMessageCalcDiatomicTwoElecsTwoCores1stDerivativesNullMatrix
106+ = "Error in am1::Am1::CalcDiatomicTwoElecsTwoCores1stDerivatives: The two elec two core diatomic matrix is NULL.\n";
107+ this->errorMessageCalcDiatomicTwoElecsTwoCores2ndDerivativesNullMatrix
108+ = "Error in am1::Am1::CalcDiatomicTwoElecsTwoCores2ndDerivatives: The two elec two core diatomic matrix is NULL.\n";
109109 this->errorMessageGetElectronicEnergyEnergyNotCalculated
110110 = "Error in am1::Am1::GetElectronicEnergy: Set electronic state is not calculated by CIS.\n";
111111 this->errorMessageGetElectronicEnergyNULLCISEnergy
--- a/src/am1/Am1D.cpp
+++ b/src/am1/Am1D.cpp
@@ -94,20 +94,20 @@ void Am1D::SetMessages(){
9494 = "Error in am1::Am1D::GetNddoRepulsionIntegral1stDerivative: Bad orbital is set.\n";
9595 this->errorMessageGetNddoRepulsionIntegral2ndDerivative
9696 = "Error in am1::Am1D::GetNddoRepulsionIntegral2ndDerivative: Bad orbital is set.\n";
97- this->errorMessageCalcTwoElecTwoCoreNullMatrix
98- = "Error in am1::Am1D::CalcTwoElecTwoCore: The two elec two core matrix is NULL.\n";
99- this->errorMessageCalcDiatomicTwoElecTwoCoreSameAtoms
100- = "Error in am1::Am1D::CalcDiatomicTwoElecTwoCore: Atom A and B is same.\n";
101- this->errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesSameAtoms
102- = "Error in am1::Am1D::CalcDiatomicTwoElecTwoCore1stDerivatives: Atom A and B is same.\n";
103- this->errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesSameAtoms
104- = "Error in am1::Am1D::CalcDiatomicTwoElecTwoCore2ndDerivatives: Atom A and B is same.\n";
105- this->errorMessageCalcDiatomicTwoElecTwoCoreNullMatrix
106- = "Error in am1::Am1D::CalcDiatomicTwoElecTwoCore: The two elec two core diatomic matrix is NULL.\n";
107- this->errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesNullMatrix
108- = "Error in am1::Am1D::CalcDiatomicTwoElecTwoCore1stDerivatives: The two elec two core diatomic matrix is NULL.\n";
109- this->errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesNullMatrix
110- = "Error in am1::Am1D::CalcDiatomicTwoElecTwoCore2ndDerivatives: The two elec two core diatomic matrix is NULL.\n";
97+ this->errorMessageCalcTwoElecsTwoCoresNullMatrix
98+ = "Error in am1::Am1D::CalcTwoElecsTwoCores: The two elec two core matrix is NULL.\n";
99+ this->errorMessageCalcDiatomicTwoElecsTwoCoresSameAtoms
100+ = "Error in am1::Am1D::CalcDiatomicTwoElecsTwoCores: Atom A and B is same.\n";
101+ this->errorMessageCalcDiatomicTwoElecsTwoCores1stDerivativesSameAtoms
102+ = "Error in am1::Am1D::CalcDiatomicTwoElecsTwoCores1stDerivatives: Atom A and B is same.\n";
103+ this->errorMessageCalcDiatomicTwoElecsTwoCores2ndDerivativesSameAtoms
104+ = "Error in am1::Am1D::CalcDiatomicTwoElecsTwoCores2ndDerivatives: Atom A and B is same.\n";
105+ this->errorMessageCalcDiatomicTwoElecsTwoCoresNullMatrix
106+ = "Error in am1::Am1D::CalcDiatomicTwoElecsTwoCores: The two elec two core diatomic matrix is NULL.\n";
107+ this->errorMessageCalcDiatomicTwoElecsTwoCores1stDerivativesNullMatrix
108+ = "Error in am1::Am1D::CalcDiatomicTwoElecsTwoCores1stDerivatives: The two elec two core diatomic matrix is NULL.\n";
109+ this->errorMessageCalcDiatomicTwoElecsTwoCores2ndDerivativesNullMatrix
110+ = "Error in am1::Am1D::CalcDiatomicTwoElecsTwoCores2ndDerivatives: The two elec two core diatomic matrix is NULL.\n";
111111 this->errorMessageGetElectronicEnergyEnergyNotCalculated
112112 = "Error in am1::Am1D::GetElectronicEnergy: Set electronic state is not calculated by CIS.\n";
113113 this->errorMessageGetElectronicEnergyNULLCISEnergy
--- a/src/base/Enums.h
+++ b/src/base/Enums.h
@@ -125,6 +125,7 @@ RENUMSTR_BEGIN( AtomType, AtomTypeStr )
125125 RENUMSTR( S, "S" )
126126 RENUMSTR( Cl, "Cl" )
127127 RENUMSTR( Ar, "Ar" )
128+ RENUMSTR( EPC, "Environmental Point Charge" )
128129 RENUMSTR( AtomType_end, "AtomType_end" )
129130 RENUMSTR_END()
130131
--- a/src/base/InputParser.cpp
+++ b/src/base/InputParser.cpp
@@ -82,6 +82,8 @@ void InputParser::SetMessages(){
8282 = "Error in base::InputParser::GetInputTerms: Input file is empty.\n";
8383 this->errorMessageNotFoundInputFile
8484 = "Error in base::InputParser::StoreInputTermsFromFile: Not found.\n";
85+ this->errorMessageNonValidTheoriesEpc
86+ = "Error in base::InputParser::ValidateEpcConditions: Theory you set is not supported for EPC. The supported theories are MNDO-sefies(MNDO, AM1, AM1D, PM3, PM3D, PDDG/PM3) only.\n";
8587 this->errorMessageNonValidTheoriesMD
8688 = "Error in base::InputParser::ValidateMdConditions: Theory you set is not supported for MD.\n";
8789 this->errorMessageNonValidExcitedStatesMD
@@ -248,6 +250,10 @@ void InputParser::SetMessages(){
248250 this->stringGeometry = "geometry";
249251 this->stringGeometryEnd = "geometry_end";
250252
253+ // Environmental Point Charge
254+ this->stringEpc = "epc";
255+ this->stringEpcEnd = "epc_end";
256+
251257 // SCF
252258 this->stringScf = "scf";
253259 this->stringScfEnd = "scf_end";
@@ -488,6 +494,22 @@ int InputParser::ParseMolecularGeometry(Molecule* molecule, vector<string>* inpu
488494 return parseIndex;
489495 }
490496
497+int InputParser::ParseEpcsConfiguration(Molecule* molecule, vector<string>* inputTerms, int parseIndex) const{
498+ parseIndex++;
499+ while((*inputTerms)[parseIndex].compare(this->stringEpcEnd) != 0){
500+ double x = atof((*inputTerms)[parseIndex+0].c_str()) * Parameters::GetInstance()->GetAngstrom2AU();
501+ double y = atof((*inputTerms)[parseIndex+1].c_str()) * Parameters::GetInstance()->GetAngstrom2AU();
502+ double z = atof((*inputTerms)[parseIndex+2].c_str()) * Parameters::GetInstance()->GetAngstrom2AU();
503+ double charge = atof((*inputTerms)[parseIndex+3].c_str());
504+ AtomType atomType = EPC;
505+ int index = molecule->GetNumberEpcs();
506+ Atom* atom = AtomFactory::Create(atomType, index, x, y, z, charge);
507+ molecule->AddEpc(atom);
508+ parseIndex += 4;
509+ }
510+ return parseIndex;
511+}
512+
491513 int InputParser::ParseConditionsSCF(vector<string>* inputTerms, int parseIndex) const{
492514 parseIndex++;
493515 while((*inputTerms)[parseIndex].compare(this->stringScfEnd) != 0){
@@ -1169,6 +1191,11 @@ void InputParser::Parse(Molecule* molecule, int argc, char *argv[]) const{
11691191 i = this->ParseMolecularGeometry(molecule, &inputTerms, i);
11701192 }
11711193
1194+ // Environmental Point Charges Configuration
1195+ if(inputTerms[i].compare(this->stringEpc) == 0){
1196+ i = this->ParseEpcsConfiguration(molecule, &inputTerms, i);
1197+ }
1198+
11721199 // scf condition
11731200 if(inputTerms[i].compare(this->stringScf) == 0){
11741201 i = this->ParseConditionsSCF(&inputTerms, i);
@@ -1249,6 +1276,7 @@ void InputParser::Parse(Molecule* molecule, int argc, char *argv[]) const{
12491276 // calculate basics and validate conditions
12501277 this->CalcMolecularBasics(molecule);
12511278 this->ValidateVdWConditions();
1279+ this->ValidateEpcConditions(*molecule);
12521280 if(Parameters::GetInstance()->RequiresCIS()){
12531281 this->ValidateCisConditions(*molecule);
12541282 }
@@ -1326,6 +1354,25 @@ void InputParser::ValidateVdWConditions() const{
13261354 }
13271355 }
13281356
1357+void InputParser::ValidateEpcConditions(const Molecule& molecule) const{
1358+ if(molecule.GetNumberEpcs()<=0){return;}
1359+ TheoryType theory = Parameters::GetInstance()->GetCurrentTheory();
1360+ // Validate theory
1361+ if(theory == MNDO ||
1362+ theory == AM1 ||
1363+ theory == AM1D ||
1364+ theory == PM3 ||
1365+ theory == PM3D ||
1366+ theory == PM3PDDG ){
1367+ }
1368+ else{
1369+ stringstream ss;
1370+ ss << this->errorMessageNonValidTheoriesEpc;
1371+ ss << this->errorMessageTheory << TheoryTypeStr(theory) << endl;
1372+ throw MolDSException(ss.str());
1373+ }
1374+}
1375+
13291376 void InputParser::ValidateCisConditions(const Molecule& molecule) const{
13301377
13311378 // direct CIS
@@ -1550,6 +1597,7 @@ void InputParser::OutputMolecularBasics(Molecule* molecule) const{
15501597 molecule->OutputConfiguration();
15511598 molecule->OutputXyzCOM();
15521599 molecule->OutputXyzCOC();
1600+ molecule->OutputEpcs();
15531601 }
15541602
15551603 void InputParser::OutputScfConditions() const{
--- a/src/base/InputParser.h
+++ b/src/base/InputParser.h
@@ -35,6 +35,7 @@ private:
3535 void SetMessages();
3636 std::string errorMessageInputFileEmpty;
3737 std::string errorMessageNotFoundInputFile;
38+ std::string errorMessageNonValidTheoriesEpc;
3839 std::string errorMessageNonValidTheoriesMD;
3940 std::string errorMessageNonValidExcitedStatesMD;
4041 std::string errorMessageNonValidExcitedStatesMC;
@@ -173,6 +174,9 @@ private:
173174 // geometry
174175 std::string stringGeometry;
175176 std::string stringGeometryEnd;
177+ // EPC
178+ std::string stringEpc;
179+ std::string stringEpcEnd;
176180 // SCF
177181 std::string stringScf;
178182 std::string stringScfEnd;
@@ -297,6 +301,7 @@ private:
297301 std::string stringFrequenciesElecState;
298302 void CalcMolecularBasics(Molecule* molecule) const;
299303 void ValidateVdWConditions() const;
304+ void ValidateEpcConditions(const Molecule& molecule) const;
300305 void ValidateCisConditions(const Molecule& molecule) const;
301306 void ValidateMdConditions(const Molecule& molecule) const;
302307 void ValidateMcConditions(const Molecule& molecule) const;
@@ -324,6 +329,7 @@ private:
324329 void StoreInputTermsFromFile(std::vector<std::string>& inputTerms, char* fileName) const;
325330 void AddInputTermsFromString(std::vector<std::string>& inputTerms, std::string str) const;
326331 int ParseMolecularGeometry(Molecule* molecule, std::vector<std::string>* inputTerms, int parseIndex) const;
332+ int ParseEpcsConfiguration(Molecule* molecule, std::vector<std::string>* inputTerms, int parseIndex) const;
327333 int ParseTheory(std::vector<std::string>* inputTerms, int parseIndex) const;
328334 int ParseConditionsSCF(std::vector<std::string>* inputTerms, int parseIndex) const;
329335 int ParseConditionsPrincipalAxes(std::vector<std::string>* inputTerms, int parseIndex) const;
--- a/src/base/MallocerFreer.h
+++ b/src/base/MallocerFreer.h
@@ -30,6 +30,7 @@ public:
3030 //1d
3131 template<typename T> void Malloc(T** matrix, size_t size1) const{
3232 if(*matrix!=NULL) return;
33+ if(size1<=0) return;
3334 double requiredMalloc = this->GetMemoryAmount<T>(size1);
3435 this->CheckLimitHeap(requiredMalloc);
3536 *matrix = new T[size1]();
@@ -58,6 +59,7 @@ public:
5859 // 2d
5960 template<typename T> void Malloc(T*** matrix, size_t size1, size_t size2) const{
6061 if(*matrix!=NULL) return;
62+ if(size1*size2<=0) return;
6163 double requiredMalloc = this->GetMemoryAmount<T*>(size1);
6264 this->CheckLimitHeap(requiredMalloc);
6365
@@ -100,6 +102,7 @@ public:
100102 // 3d
101103 template<typename T> void Malloc(T**** matrix, size_t size1, size_t size2, size_t size3) const{
102104 if(*matrix!=NULL) return;
105+ if(size1*size2*size3<=0) return;
103106 double requiredMalloc = this->GetMemoryAmount<T**>(size1);
104107 this->CheckLimitHeap(requiredMalloc);
105108
@@ -144,6 +147,7 @@ public:
144147 //4d
145148 template<typename T> void Malloc(T***** matrix, size_t size1, size_t size2, size_t size3, size_t size4) const{
146149 if(*matrix!=NULL) return;
150+ if(size1*size2*size3*size4<=0) return;
147151 double requiredMalloc = this->GetMemoryAmount<T***>(size1);
148152 this->CheckLimitHeap(requiredMalloc);
149153
@@ -189,6 +193,7 @@ public:
189193 //5d
190194 template<typename T> void Malloc(T****** matrix, size_t size1, size_t size2, size_t size3, size_t size4, size_t size5) const{
191195 if(*matrix!=NULL) return;
196+ if(size1*size2*size3*size4*size5<=0) return;
192197 double requiredMalloc = this->GetMemoryAmount<T****>(size1);
193198 this->CheckLimitHeap(requiredMalloc);
194199
@@ -235,6 +240,7 @@ public:
235240 //6d
236241 template<typename T> void Malloc(T******* matrix, size_t size1, size_t size2, size_t size3, size_t size4, size_t size5, size_t size6) const{
237242 if(*matrix!=NULL) return;
243+ if(size1*size2*size3*size4*size5*size6<=0) return;
238244 double requiredMalloc = this->GetMemoryAmount<T*****>(size1);
239245 this->CheckLimitHeap(requiredMalloc);
240246
@@ -282,6 +288,7 @@ public:
282288 //7d
283289 template<typename T> void Malloc(T******** matrix, size_t size1, size_t size2, size_t size3, size_t size4, size_t size5, size_t size6, size_t size7) const{
284290 if(*matrix!=NULL) return;
291+ if(size1*size2*size3*size4*size5*size6*size7<=0) return;
285292 double requiredMalloc = this->GetMemoryAmount<T******>(size1);
286293 this->CheckLimitHeap(requiredMalloc);
287294
--- a/src/base/Molecule.cpp
+++ b/src/base/Molecule.cpp
@@ -57,15 +57,30 @@ Molecule::Molecule(const Molecule& rhs){
5757 Molecule& Molecule::operator=(const Molecule& rhs){
5858 double* oldXyzCOM = this->xyzCOM;
5959 double* oldXyzCOC = this->xyzCOC;
60- double** oldDistanceMatrix = this->distanceMatrix;
60+ double** oldDistanceAtoms = this->distanceAtoms;
61+ double** oldDistanceEpcs = this->distanceEpcs;
62+ double** oldDistanceAtomsEpcs = this->distanceAtomsEpcs;
6163 vector<Atom*>* oldAtomVect = this->atomVect;
64+ vector<Atom*>* oldEpcVect = this->epcVect;
6265 this->CopyInitialize(rhs);
63- this->Finalize(&oldAtomVect, &oldXyzCOM, &oldXyzCOC, &oldDistanceMatrix);
66+ this->Finalize(&oldAtomVect,
67+ &oldEpcVect,
68+ &oldXyzCOM,
69+ &oldXyzCOC,
70+ &oldDistanceAtoms,
71+ &oldDistanceEpcs,
72+ &oldDistanceAtomsEpcs);
6473 return *this;
6574 }
6675
6776 Molecule::~Molecule(){
68- this->Finalize(&this->atomVect, &this->xyzCOM, &this->xyzCOC, &this->distanceMatrix);
77+ this->Finalize(&this->atomVect,
78+ &this->epcVect,
79+ &this->xyzCOM,
80+ &this->xyzCOC,
81+ &this->distanceAtoms,
82+ &this->distanceEpcs,
83+ &this->distanceAtomsEpcs);
6984 //this->OutputLog("molecule deleted\n");
7085 }
7186
@@ -78,48 +93,101 @@ void Molecule::CopyInitialize(const Molecule& rhs){
7893 this->totalNumberAOs = rhs.totalNumberAOs;
7994 this->totalNumberValenceElectrons = rhs.totalNumberValenceElectrons;
8095 this->totalCoreMass = rhs.totalCoreMass;
81- for(int i=0; i<rhs.atomVect->size(); i++){
82- Atom* atom = (*rhs.atomVect)[i];
83- this->atomVect->push_back(AtomFactory::Create(atom->GetAtomType(),
84- atom->GetIndex(),
85- atom->GetXyz()[XAxis],
86- atom->GetXyz()[YAxis],
87- atom->GetXyz()[ZAxis],
88- atom->GetPxyz()[XAxis],
89- atom->GetPxyz()[YAxis],
90- atom->GetPxyz()[ZAxis]));
91- (*this->atomVect)[i]->SetFirstAOIndex(atom->GetFirstAOIndex());
92- }
93- int atomNum = this->atomVect->size();
94- MallocerFreer::GetInstance()->Malloc<double>(&this->distanceMatrix, atomNum, atomNum);
95- for(int i=0; i<atomNum; i++){
96- for(int j=0; j<atomNum; j++){
97- this->distanceMatrix[i][j] = rhs.distanceMatrix[i][j];
96+ if(rhs.atomVect != NULL){
97+ int atomNum = rhs.atomVect->size();
98+ for(int i=0; i<atomNum; i++){
99+ Atom* atom = (*rhs.atomVect)[i];
100+ this->atomVect->push_back(AtomFactory::Create(atom->GetAtomType(),
101+ atom->GetIndex(),
102+ atom->GetXyz()[XAxis],
103+ atom->GetXyz()[YAxis],
104+ atom->GetXyz()[ZAxis],
105+ atom->GetPxyz()[XAxis],
106+ atom->GetPxyz()[YAxis],
107+ atom->GetPxyz()[ZAxis]));
108+ (*this->atomVect)[i]->SetFirstAOIndex(atom->GetFirstAOIndex());
109+ }
110+ MallocerFreer::GetInstance()->Malloc<double>(&this->distanceAtoms, atomNum, atomNum);
111+ for(int i=0; i<atomNum; i++){
112+ for(int j=0; j<atomNum; j++){
113+ this->distanceAtoms[i][j] = rhs.distanceAtoms[i][j];
114+ }
115+ }
116+ }
117+ if(rhs.epcVect != NULL){
118+ int epcNum = rhs.epcVect->size();
119+ for(int i=0; i<epcNum; i++){
120+ Atom* epc = (*rhs.epcVect)[i];
121+ this->epcVect->push_back(
122+ AtomFactory::Create(EPC,
123+ epc->GetIndex(),
124+ epc->GetXyz()[XAxis],
125+ epc->GetXyz()[YAxis],
126+ epc->GetXyz()[ZAxis],
127+ epc->GetPxyz()[XAxis],
128+ epc->GetPxyz()[YAxis],
129+ epc->GetPxyz()[ZAxis],
130+ epc->GetCoreCharge()));
131+ (*this->epcVect)[i]->SetFirstAOIndex(epc->GetFirstAOIndex());
132+ }
133+ MallocerFreer::GetInstance()->Malloc<double>(&this->distanceEpcs, epcNum, epcNum);
134+ for(int i=0; i<epcNum; i++){
135+ for(int j=0; j<epcNum; j++){
136+ this->distanceEpcs[i][j] = rhs.distanceEpcs[i][j];
137+ }
138+ }
139+ }
140+ if(rhs.atomVect != NULL && rhs.epcVect != NULL){
141+ int atomNum = rhs.atomVect->size();
142+ int epcNum = rhs.epcVect->size();
143+ MallocerFreer::GetInstance()->Malloc<double>(&this->distanceAtomsEpcs, atomNum, epcNum);
144+ for(int i=0; i<atomNum; i++){
145+ for(int j=0; j<epcNum; j++){
146+ this->distanceAtomsEpcs[i][j] = rhs.distanceAtomsEpcs[i][j];
147+ }
98148 }
99149 }
100-
101150 }
102151
103152 void Molecule::Initialize(){
104153 this->SetMessages();
105- this->xyzCOM = NULL;
106- this->xyzCOC = NULL;
107- this->distanceMatrix = NULL;
108- this->atomVect = NULL;
154+ this->xyzCOM = NULL;
155+ this->xyzCOC = NULL;
156+ this->distanceAtoms = NULL;
157+ this->distanceEpcs = NULL;
158+ this->distanceAtomsEpcs = NULL;
159+ this->atomVect = NULL;
160+ this->epcVect = NULL;
109161 try{
110162 this->atomVect = new vector<Atom*>;
163+ this->epcVect = new vector<Atom*>;
111164 MallocerFreer::GetInstance()->Malloc<double>(&this->xyzCOM, CartesianType_end);
112165 MallocerFreer::GetInstance()->Malloc<double>(&this->xyzCOC, CartesianType_end);
113166 }
114167 catch(exception ex){
115- this->Finalize(&this->atomVect, &this->xyzCOM, &this->xyzCOC, &this->distanceMatrix);
168+ this->Finalize(&this->atomVect,
169+ &this->epcVect,
170+ &this->xyzCOM,
171+ &this->xyzCOC,
172+ &this->distanceAtoms,
173+ &this->distanceEpcs,
174+ &this->distanceAtomsEpcs);
116175 throw MolDSException(ex.what());
117176 }
118177 }
119178
120-void Molecule::Finalize(vector<Atom*>** atomVect, double** xyzCOM, double** xyzCOC, double*** distanceMatrix){
121- int atomNum = (*atomVect)->size();
179+void Molecule::Finalize(vector<Atom*>** atomVect,
180+ vector<Atom*>** epcVect,
181+ double** xyzCOM,
182+ double** xyzCOC,
183+ double*** distanceAtoms,
184+ double*** distanceEpcs,
185+ double*** distanceAtomsEpcs){
186+ MallocerFreer::GetInstance()->Free<double>(xyzCOM, CartesianType_end);
187+ MallocerFreer::GetInstance()->Free<double>(xyzCOC, CartesianType_end);
188+ int atomNum=0;
122189 if(*atomVect != NULL){
190+ atomNum = (*atomVect)->size();
123191 for(int i=0; i<atomNum; i++){
124192 if((**atomVect)[i] != NULL){
125193 delete (**atomVect)[i];
@@ -130,16 +198,37 @@ void Molecule::Finalize(vector<Atom*>** atomVect, double** xyzCOM, double** xyzC
130198 delete *atomVect;
131199 *atomVect = NULL;
132200 //this->OutputLog("atomVect deleted\n");
201+ MallocerFreer::GetInstance()->Free<double>(distanceAtoms, atomNum, atomNum);
202+ }
203+ int epcNum = 0;
204+ if(*epcVect != NULL){
205+ epcNum = (*epcVect)->size();
206+ for(int i=0; i<epcNum; i++){
207+ if((**epcVect)[i] != NULL){
208+ delete (**epcVect)[i];
209+ (**epcVect)[i] = NULL;
210+ }
211+ }
212+ (*epcVect)->clear();
213+ delete *epcVect;
214+ *epcVect = NULL;
215+ //this->OutputLog("epcVect deleted\n");
216+ MallocerFreer::GetInstance()->Free<double>(distanceEpcs, epcNum, epcNum);
217+ }
218+ if(*atomVect != NULL && *epcVect != NULL){
219+ atomNum = (*atomVect)->size();
220+ epcNum = (*epcVect)->size();
221+ MallocerFreer::GetInstance()->Free<double>(distanceAtomsEpcs, atomNum, epcNum);
133222 }
134- MallocerFreer::GetInstance()->Free<double>(xyzCOM, CartesianType_end);
135- MallocerFreer::GetInstance()->Free<double>(xyzCOC, CartesianType_end);
136- MallocerFreer::GetInstance()->Free<double>(distanceMatrix, atomNum, atomNum);
137223 }
138224
139225 void Molecule::SetMessages(){
140226 this->errorMessageGetAtomNull = "Error in base::Molecule::GetAtom: atomVect is NULL.\n";
227+ this->errorMessageGetEPCNull = "Error in base::Molecule::GetEnviromentalPointCharge: enviromentalPointChargeVect is NULL.\n";
141228 this->errorMessageAddAtomNull = "Error in base::Molecule::AddAtom: atomVect is NULL.\n";
229+ this->errorMessageAddEPCNull = "Error in base::Molecule::AddEnviromentalPointCharge: enviromentalPointChargeVect is NULL.\n";
142230 this->errorMessageGetNumberAtomsNull = "Error in base::Molecule::GetNumberAtoms: atomVect is NULL.\n";
231+ this->errorMessageGetNumberEPCsNull = "Error in base::Molecule::GetNumberEnviromentalPointChargess: epcVect is NULL.\n";
143232 this->errorMessageGetXyzCOMNull = "Error in base::Molecule::GetXyzCOM: xyzCOM is NULL.\n";
144233 this->errorMessageGetXyzCOCNull = "Error in base::Molecule::GetXyzCOC: xyzCOC is NULL.\n";
145234 this->errorMessageCalcXyzCOMNull = "Error in base::Molecule::CalcXyzCOM: xyzCOM is NULL.\n";
@@ -151,6 +240,9 @@ void Molecule::SetMessages(){
151240 this->messageAtomCoordinates = "\tAtom coordinates:";
152241 this->messageAtomMomenta = "\tAtom momenta:";
153242 this->messageAtomMomentaTitle = "\t\t\t| i-th | atom type | px[a.u.] | py[a.u.] | pz[a.u.] |\t\t| px[u] | py[u] | pz[u] | [u] = [(g/Mol)*(angst/fs)]\n";
243+ this->messageEpcConfiguration = "\tEnvironmental Point Charge(EPC) configuration:\n";
244+ this->messageEpcCoordinates = "\tEPC coordinates:";
245+ this->messageEpcCoordinatesTitle = "\t\t\t\t| i-th | charge[a.u.] | x[a.u.] | y[a.u.] | z[a.u.] |\t\t| x[angst.] | y[angst.] | z[angst.] |\n";
154246 this->messageCOM = "\tCenter of Mass:";
155247 this->messageCOC = "\tCenter of Core:";
156248 this->messageCOMTitle = "\t\t\t| x[a.u.] | y[a.u.] | z[a.u.] |\t\t| x[angst.] | y[angst.] | z[angst.] |\n";
@@ -182,7 +274,13 @@ void Molecule::AddAtom(Atom* atom){
182274 if(this->atomVect==NULL) throw MolDSException(this->errorMessageAddAtomNull);
183275 #endif
184276 this->atomVect->push_back(atom);
185-
277+}
278+
279+void Molecule::AddEpc(Atom* epc){
280+#ifdef MOLDS_DBG
281+ if(this->epcVect==NULL) throw MolDSException(this->errorMessageAddEPCNull);
282+#endif
283+ this->epcVect->push_back(epc);
186284 }
187285
188286 double const* Molecule::GetXyzCOM() const{
@@ -251,9 +349,9 @@ void Molecule::CalcXyzCOC(){
251349 }
252350 }
253351
254-void Molecule::CalcDistanceMatrix(){
255- if(this->distanceMatrix==NULL){
256- MallocerFreer::GetInstance()->Malloc<double>(&this->distanceMatrix, this->atomVect->size(), this->atomVect->size());
352+void Molecule::CalcDistanceAtoms(){
353+ if(this->distanceAtoms==NULL){
354+ MallocerFreer::GetInstance()->Malloc<double>(&this->distanceAtoms, this->atomVect->size(), this->atomVect->size());
257355 }
258356 for(int a=0; a<this->atomVect->size(); a++){
259357 const Atom& atomA = *(*this->atomVect)[a];
@@ -263,8 +361,45 @@ void Molecule::CalcDistanceMatrix(){
263361 distance = sqrt( pow(atomA.GetXyz()[0] - atomB.GetXyz()[0], 2.0)
264362 +pow(atomA.GetXyz()[1] - atomB.GetXyz()[1], 2.0)
265363 +pow(atomA.GetXyz()[2] - atomB.GetXyz()[2], 2.0) );
266- this->distanceMatrix[a][b] = distance;
267- this->distanceMatrix[b][a] = distance;
364+ this->distanceAtoms[a][b] = distance;
365+ this->distanceAtoms[b][a] = distance;
366+ }
367+ }
368+}
369+
370+void Molecule::CalcDistanceEpcs(){
371+ if(this->epcVect == NULL){return;}
372+ if(this->distanceEpcs==NULL){
373+ MallocerFreer::GetInstance()->Malloc<double>(&this->distanceEpcs, this->epcVect->size(), this->epcVect->size());
374+ }
375+ for(int a=0; a<this->epcVect->size(); a++){
376+ const Atom& epcA = *(*this->epcVect)[a];
377+ for(int b=a; b<this->epcVect->size(); b++){
378+ const Atom& epcB = *(*this->epcVect)[b];
379+ double distance=0.0;
380+ distance = sqrt( pow(epcA.GetXyz()[0] - epcB.GetXyz()[0], 2.0)
381+ +pow(epcA.GetXyz()[1] - epcB.GetXyz()[1], 2.0)
382+ +pow(epcA.GetXyz()[2] - epcB.GetXyz()[2], 2.0) );
383+ this->distanceEpcs[a][b] = distance;
384+ this->distanceEpcs[b][a] = distance;
385+ }
386+ }
387+}
388+
389+void Molecule::CalcDistanceAtomsEpcs(){
390+ if(this->epcVect == NULL){return;}
391+ if(this->distanceAtomsEpcs==NULL){
392+ MallocerFreer::GetInstance()->Malloc<double>(&this->distanceAtomsEpcs, this->atomVect->size(), this->epcVect->size());
393+ }
394+ for(int a=0; a<this->atomVect->size(); a++){
395+ const Atom& atom = *(*this->atomVect)[a];
396+ for(int b=0; b<this->epcVect->size(); b++){
397+ const Atom& epc = *(*this->epcVect)[b];
398+ double distance=0.0;
399+ distance = sqrt( pow(atom.GetXyz()[0] - epc.GetXyz()[0], 2.0)
400+ +pow(atom.GetXyz()[1] - epc.GetXyz()[1], 2.0)
401+ +pow(atom.GetXyz()[2] - epc.GetXyz()[2], 2.0) );
402+ this->distanceAtomsEpcs[a][b] = distance;
268403 }
269404 }
270405 }
@@ -279,7 +414,9 @@ void Molecule::CalcBasics(){
279414 void Molecule::CalcBasicsConfiguration(){
280415 this->CalcXyzCOM();
281416 this->CalcXyzCOC();
282- this->CalcDistanceMatrix();
417+ this->CalcDistanceAtoms();
418+ this->CalcDistanceEpcs();
419+ this->CalcDistanceAtomsEpcs();
283420 }
284421
285422 void Molecule::CalcTotalNumberAOs(){
@@ -347,6 +484,27 @@ void Molecule::OutputMomenta() const{
347484 this->OutputLog("\n");
348485 }
349486
487+void Molecule::OutputEpcs() const{
488+ if(this->epcVect == NULL || this->epcVect->size() <= 0) {return;}
489+ double ang2AU = Parameters::GetInstance()->GetAngstrom2AU();
490+ this->OutputLog(this->messageEpcConfiguration);
491+ this->OutputLog(this->messageEpcCoordinatesTitle);
492+ for(int a=0; a<this->epcVect->size(); a++){
493+ const Atom& atom = *(*this->epcVect)[a];
494+ this->OutputLog(boost::format("%s\t%d\t%e\t%e\t%e\t%e\t\t%e\t%e\t%e\n")
495+ % this->messageEpcCoordinates
496+ % a
497+ % atom.GetCoreCharge()
498+ % atom.GetXyz()[0]
499+ % atom.GetXyz()[1]
500+ % atom.GetXyz()[2]
501+ % (atom.GetXyz()[0]/ang2AU)
502+ % (atom.GetXyz()[1]/ang2AU)
503+ % (atom.GetXyz()[2]/ang2AU));
504+ }
505+ this->OutputLog("\n");
506+}
507+
350508 void Molecule::OutputXyzCOM() const{
351509 double ang2AU = Parameters::GetInstance()->GetAngstrom2AU();
352510 this->OutputLog(this->messageCOMTitle);
@@ -673,14 +831,6 @@ void Molecule::OutputTranslatingConditions(double const* translatingDifference)
673831 % (translatingDifference[2]/angst2AU));
674832 }
675833
676-//double Molecule::GetDistanceAtoms(int indexAtomA, int indexAtomB) const{
677-// return this->distanceMatrix[indexAtomA][indexAtomB];
678-//}
679-//
680-//double Molecule::GetDistanceAtoms(const Atom& atomA, const Atom& atomB) const{
681-// return this->GetDistanceAtoms(atomA.GetIndex(), atomB.GetIndex());
682-//}
683-
684834 void Molecule::SynchronizeConfigurationTo(const Molecule& ref){
685835 for(int a=0; a<this->GetNumberAtoms(); a++){
686836 Atom& atom = *this->GetAtom(a);
--- a/src/base/Molecule.h
+++ b/src/base/Molecule.h
@@ -38,7 +38,20 @@ public:
3838 #endif
3939 return (*this->atomVect)[atomIndex];
4040 }
41+ inline int GetNumberEpcs() const{
42+#ifdef MOLDS_DBG
43+ if(this->epcVect==NULL) throw MolDS_base::MolDSException(this->errorMessageGetNumberEPCsNull);
44+#endif
45+ return this->epcVect->size();
46+ }
47+ inline MolDS_base_atoms::Atom* GetEpc(int epcIndex) const{
48+#ifdef MOLDS_DBG
49+ if(this->epcVect==NULL) throw MolDS_base::MolDSException(this->errorMessageGetEPCNull);
50+#endif
51+ return (*this->epcVect)[epcIndex];
52+ }
4153 void AddAtom(MolDS_base_atoms::Atom* atom);
54+ void AddEpc(MolDS_base_atoms::Atom* epc);
4255 double const* GetXyzCOM() const;
4356 double const* GetXyzCOC() const;
4457 void CalcBasics();
@@ -51,12 +64,19 @@ public:
5164 void OutputTotalNumberAtomsAOsValenceelectrons() const;
5265 void OutputConfiguration() const;
5366 void OutputMomenta() const;
67+ void OutputEpcs() const;
5468 void CalcPrincipalAxes();
5569 void Rotate();
5670 void Translate();
57- double GetDistanceAtoms(int indexAtomA, int indexAtomB) const{return this->distanceMatrix[indexAtomA][indexAtomB];};
71+ double GetDistanceAtoms(int indexAtomA, int indexAtomB) const{return this->distanceAtoms[indexAtomA][indexAtomB];};
5872 double GetDistanceAtoms(const MolDS_base_atoms::Atom& atomA,
5973 const MolDS_base_atoms::Atom& atomB) const{return this->GetDistanceAtoms(atomA.GetIndex(), atomB.GetIndex());};
74+ double GetDistanceEpcs(int indexEpcA, int indexEpcB) const{return this->distanceEpcs[indexEpcA][indexEpcB];};
75+ double GetDistanceEpcs(const MolDS_base_atoms::Atom& epcA,
76+ const MolDS_base_atoms::Atom& epcB) const{return this->GetDistanceEpcs(epcA.GetIndex(), epcB.GetIndex());};
77+ double GetDistanceAtomEpc(int indexAtom, int indexEpc) const{return this->distanceAtomsEpcs[indexAtom][indexEpc];};
78+ double GetDistanceAtomEpc(const MolDS_base_atoms::Atom& atom,
79+ const MolDS_base_atoms::Atom& epc) const{return this->GetDistanceAtomEpc(atom.GetIndex(), epc.GetIndex());};
6080 void SynchronizeConfigurationTo (const Molecule& ref);
6181 void SynchronizeMomentaTo (const Molecule& ref);
6282 void SynchronizePhaseSpacePointTo(const Molecule& ref);
@@ -65,22 +85,33 @@ public:
6585 void BroadcastPhaseSpacePointToAllProcesses(int root) const;
6686 private:
6787 std::vector<MolDS_base_atoms::Atom*>* atomVect;
88+ std::vector<MolDS_base_atoms::Atom*>* epcVect; // Vector of Environmental Point Charges
6889 double* xyzCOM; // x, y, z coordinates of Center of Mass;
6990 double* xyzCOC; // x, y, z coordinates of Center of Core;
70- double** distanceMatrix; // distance between each atom;
91+ double** distanceAtoms; // distance between each atom;
92+ double** distanceEpcs; // distance between each environmental point charge;
93+ double** distanceAtomsEpcs;// distance between each atom and environmental point charge;
7194 int totalNumberAOs;
7295 int totalNumberValenceElectrons;
7396 double totalCoreMass;
7497 void Initialize();
7598 void CopyInitialize(const Molecule& rhs);
76- void Finalize(std::vector<MolDS_base_atoms::Atom*>** atomVect, double** xyzCOM, double** xyzCOC, double*** distanceMatrix);
99+ void Finalize(std::vector<MolDS_base_atoms::Atom*>** atomVect,
100+ std::vector<MolDS_base_atoms::Atom*>** epcVect,
101+ double** xyzCOM,
102+ double** xyzCOC,
103+ double*** distanceAtoms,
104+ double*** distanceEpcs,
105+ double*** distanceAtomsEpcs);
77106 void SetMessages();
78107 void CalcTotalNumberValenceElectrons();
79108 void CalcTotalNumberAOs();
80109 void CalcTotalCoreMass();
81110 void CalcXyzCOM();
82111 void CalcXyzCOC();
83- void CalcDistanceMatrix();
112+ void CalcDistanceAtoms();
113+ void CalcDistanceEpcs();
114+ void CalcDistanceAtomsEpcs();
84115 void CalcInertiaTensor(double** inertiaTensor,
85116 double const* inertiaTensorOrigin);
86117 void FreeInertiaTensorMoments(double*** inertiaTensor,
@@ -98,8 +129,11 @@ private:
98129 MolDS_base::EularAngle rotatingEularAngles)const;
99130 void OutputTranslatingConditions(double const* translatingDifference) const;
100131 std::string errorMessageGetAtomNull;
132+ std::string errorMessageGetEPCNull;
101133 std::string errorMessageAddAtomNull;
134+ std::string errorMessageAddEPCNull;
102135 std::string errorMessageGetNumberAtomsNull;
136+ std::string errorMessageGetNumberEPCsNull;
103137 std::string errorMessageGetXyzCOMNull;
104138 std::string errorMessageGetXyzCOCNull;
105139 std::string errorMessageCalcXyzCOMNull;
@@ -111,6 +145,9 @@ private:
111145 std::string messageAtomCoordinatesTitle;
112146 std::string messageAtomMomenta;
113147 std::string messageAtomMomentaTitle;
148+ std::string messageEpcConfiguration;
149+ std::string messageEpcCoordinates;
150+ std::string messageEpcCoordinatesTitle;
114151 std::string messageCOM;
115152 std::string messageCOC;
116153 std::string messageCOMTitle;
--- a/src/base/atoms/Atom.h
+++ b/src/base/atoms/Atom.h
@@ -47,6 +47,7 @@ public:
4747 double GetBondingParameter(MolDS_base::TheoryType theory,
4848 MolDS_base::OrbitalType orbital) const;
4949 inline double GetCoreCharge() const{return this->coreCharge;}
50+ inline void SetCoreCharge(double charge) {this->coreCharge=charge;}
5051 inline int GetFirstAOIndex() const{return this->firstAOIndex;}
5152 inline void SetFirstAOIndex(int firstAOIndex){this->firstAOIndex = firstAOIndex;}
5253 inline int GetLastAOIndex() const{return this->firstAOIndex + this->valence.size()-1;}
--- /dev/null
+++ b/src/base/atoms/mm/EnvironmentalPointCharge.cpp
@@ -0,0 +1,52 @@
1+//************************************************************************//
2+// Copyright (C) 2011-2012 Mikiya Fujii //
3+// //
4+// This file is part of MolDS. //
5+// //
6+// MolDS is free software: you can redistribute it and/or modify //
7+// it under the terms of the GNU General Public License as published by //
8+// the Free Software Foundation, either version 3 of the License, or //
9+// (at your option) any later version. //
10+// //
11+// MolDS is distributed in the hope that it will be useful, //
12+// but WITHOUT ANY WARRANTY; without even the implied warranty of //
13+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
14+// GNU General Public License for more details. //
15+// //
16+// You should have received a copy of the GNU General Public License //
17+// along with MolDS. If not, see <http://www.gnu.org/licenses/>. //
18+//************************************************************************//
19+#include<stdio.h>
20+#include<stdlib.h>
21+#include<iostream>
22+#include<sstream>
23+#include<math.h>
24+#include<vector>
25+#include<boost/format.hpp>
26+#include"../../Enums.h"
27+#include"../../Uncopyable.h"
28+#include"../../PrintController.h"
29+#include"../../MolDSException.h"
30+#include"../../MallocerFreer.h"
31+#include"../../../mpi/MpiInt.h"
32+#include"../../../mpi/MpiProcess.h"
33+#include"../../EularAngle.h"
34+#include"../../Parameters.h"
35+#include"../../RealSphericalHarmonicsIndex.h"
36+#include"../Atom.h"
37+#include"EnvironmentalPointCharge.h"
38+using namespace std;
39+using namespace MolDS_base;
40+namespace MolDS_base_atoms_mm{
41+EnvironmentalPointCharge::EnvironmentalPointCharge(int index) : MolDS_base_atoms::Atom(index){
42+ this->SetAtomicParameters();
43+}
44+
45+void EnvironmentalPointCharge::SetAtomicParameters(){
46+ this->atomType = EPC;
47+ this->valence.push_back(s);
48+ for(int i=0; i<this->valence.size();i++){
49+ this->realSphericalHarmonicsIndeces.push_back(new RealSphericalHarmonicsIndex(this->valence[i]));
50+ }
51+}
52+}
--- /dev/null
+++ b/src/base/atoms/mm/EnvironmentalPointCharge.h
@@ -0,0 +1,30 @@
1+//************************************************************************//
2+// Copyright (C) 2011-2012 Mikiya Fujii //
3+// //
4+// This file is part of MolDS. //
5+// //
6+// MolDS is free software: you can redistribute it and/or modify //
7+// it under the terms of the GNU General Public License as published by //
8+// the Free Software Foundation, either version 3 of the License, or //
9+// (at your option) any later version. //
10+// //
11+// MolDS is distributed in the hope that it will be useful, //
12+// but WITHOUT ANY WARRANTY; without even the implied warranty of //
13+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
14+// GNU General Public License for more details. //
15+// //
16+// You should have received a copy of the GNU General Public License //
17+// along with MolDS. If not, see <http://www.gnu.org/licenses/>. //
18+//************************************************************************//
19+#ifndef INCLUDED_ENVIRONMENTAL_POINT_CHARGE
20+#define INCLUDED_ENVIRONMENTAL_POINT_CHARGE
21+namespace MolDS_base_atoms_mm{
22+class EnvironmentalPointCharge : public MolDS_base_atoms::Atom {
23+public:
24+ EnvironmentalPointCharge(int index);
25+private:
26+ EnvironmentalPointCharge();
27+ void SetAtomicParameters();
28+};
29+}
30+#endif
--- a/src/base/factories/AtomFactory.cpp
+++ b/src/base/factories/AtomFactory.cpp
@@ -38,14 +38,17 @@
3838 #include"../atoms/Natom.h"
3939 #include"../atoms/Oatom.h"
4040 #include"../atoms/Satom.h"
41+#include"../atoms/mm/EnvironmentalPointCharge.h"
4142 #include"AtomFactory.h"
4243 using namespace std;
4344 using namespace MolDS_base;
4445 using namespace MolDS_base_atoms;
46+using namespace MolDS_base_atoms_mm;
4547 namespace MolDS_base_factories{
4648
47-string AtomFactory::errorMessageNotEnableAtom = "Error in base::AtomFactory::Create: Not Enable AtomType is set.";
48-string AtomFactory::errorMessageAtomType = "\tatom type = ";
49+string AtomFactory::errorMessageNotEnableAtom = "Error in base::AtomFactory::Create: Not Enable AtomType is set.";
50+string AtomFactory::errorMessageNotEnvironmentalPointCharge = "Error in base::AtomFactory::Create: Not Environmental point charge is set.";
51+string AtomFactory::errorMessageAtomType = "\tatom type = ";
4952
5053 Atom* AtomFactory::Create(AtomType atomType, int index, double x, double y, double z, double px, double py, double pz){
5154 Atom* atom=NULL;
@@ -78,12 +81,36 @@ Atom* AtomFactory::Create(AtomType atomType, int index, double x, double y, doub
7881 return atom;
7982 }
8083
84+Atom* AtomFactory::Create(AtomType atomType, int index, double x, double y, double z, double px, double py, double pz, double charge){
85+ Atom* atom=NULL;
86+ if(atomType == EPC){
87+ atom = new EnvironmentalPointCharge(index);
88+ }
89+ else{
90+ stringstream ss;
91+ ss << AtomFactory::errorMessageNotEnvironmentalPointCharge << endl;
92+ ss << AtomFactory::errorMessageAtomType << AtomTypeStr(atomType) << endl;
93+ throw MolDSException(ss.str());
94+ }
95+ atom->SetXyz(x, y, z);
96+ atom->SetPxyz(px, py, pz);
97+ atom->SetCoreCharge(charge);
98+ return atom;
99+}
100+
81101 Atom* AtomFactory::Create(AtomType atomType, int index, double x, double y, double z){
82102 double px=0.0;
83103 double py=0.0;
84104 double pz=0.0;
85105 return AtomFactory::Create(atomType, index, x, y, z, px, py, pz);
86106 }
107+
108+Atom* AtomFactory::Create(AtomType atomType, int index, double x, double y, double z, double charge){
109+ double px=0.0;
110+ double py=0.0;
111+ double pz=0.0;
112+ return AtomFactory::Create(atomType, index, x, y, z, px, py, pz, charge);
113+}
87114 }
88115
89116
--- a/src/base/factories/AtomFactory.h
+++ b/src/base/factories/AtomFactory.h
@@ -35,11 +35,27 @@ public:
3535 int index,
3636 double x,
3737 double y,
38+ double z,
39+ double px,
40+ double py,
41+ double pz,
42+ double charge);
43+ static MolDS_base_atoms::Atom* Create(MolDS_base::AtomType atomType,
44+ int index,
45+ double x,
46+ double y,
3847 double z);
48+ static MolDS_base_atoms::Atom* Create(MolDS_base::AtomType atomType,
49+ int index,
50+ double x,
51+ double y,
52+ double z,
53+ double charge);
3954 private:
4055 AtomFactory();
4156 ~AtomFactory();
4257 static std::string errorMessageNotEnableAtom;
58+ static std::string errorMessageNotEnvironmentalPointCharge;
4359 static std::string errorMessageAtomType;
4460 };
4561
--- a/src/cndo/Cndo2.cpp
+++ b/src/cndo/Cndo2.cpp
@@ -82,7 +82,7 @@ Cndo2::Cndo2(){
8282 this->atomicElectronPopulationCIS = NULL;
8383 this->atomicUnpairedPopulationCIS = NULL;
8484 this->overlapAOs = NULL;
85- this->twoElecTwoCore = NULL;
85+ this->twoElecsTwoAtomCores = NULL;
8686 this->cartesianMatrix = NULL;
8787 this->electronicTransitionDipoleMoments = NULL;
8888 this->coreDipoleMoment = NULL;
@@ -533,7 +533,7 @@ void Cndo2::DoSCF(bool requiresGuess){
533533 this->CalcGammaAB(this->gammaAB, *this->molecule);
534534 this->CalcOverlapAOs(this->overlapAOs, *this->molecule);
535535 this->CalcCartesianMatrixByGTOExpansion(this->cartesianMatrix, *this->molecule, STO6G);
536- this->CalcTwoElecTwoCore(this->twoElecTwoCore, *this->molecule);
536+ this->CalcTwoElecsTwoCores(this->twoElecsTwoAtomCores, *this->molecule);
537537
538538 // SCF
539539 double rmsDensity=0.0;
@@ -556,7 +556,7 @@ void Cndo2::DoSCF(bool requiresGuess){
556556 this->gammaAB,
557557 this->orbitalElectronPopulation,
558558 this->atomicElectronPopulation,
559- this->twoElecTwoCore,
559+ this->twoElecsTwoAtomCores,
560560 isGuess);
561561
562562 // diagonalization of the Fock matrix
@@ -744,8 +744,8 @@ double const* const* Cndo2::GetForce(int elecState){
744744 return this->matrixForce[0];
745745 }
746746
747-void Cndo2::CalcTwoElecTwoCore(double****** twoElecTwoCore,
748- const Molecule& molecule) const{
747+void Cndo2::CalcTwoElecsTwoCores(double****** twoElecsTwoAtomCores,
748+ const Molecule& molecule) const{
749749 // do nothing for CNDO, INDO, and ZINDO/S.
750750 // two electron two core integrals are not needed for CNDO, INDO, and ZINDO/S.
751751 }
@@ -1214,7 +1214,7 @@ void Cndo2::CalcElecSCFEnergy(double* elecSCFEnergy,
12141214 this->gammaAB,
12151215 this->orbitalElectronPopulation,
12161216 this->atomicElectronPopulation,
1217- this->twoElecTwoCore,
1217+ this->twoElecsTwoAtomCores,
12181218 isGuess);
12191219 this->CalcFockMatrix(hMatrix,
12201220 molecule,
@@ -1222,7 +1222,7 @@ void Cndo2::CalcElecSCFEnergy(double* elecSCFEnergy,
12221222 this->gammaAB,
12231223 dammyOrbitalElectronPopulation,
12241224 dammyAtomicElectronPopulation,
1225- this->twoElecTwoCore,
1225+ this->twoElecsTwoAtomCores,
12261226 isGuess);
12271227
12281228 for(int i=0; i<totalNumberAOs; i++){
@@ -1386,7 +1386,7 @@ void Cndo2::CalcFockMatrix(double** fockMatrix,
13861386 double const* const* gammaAB,
13871387 double const* const* orbitalElectronPopulation,
13881388 double const* atomicElectronPopulation,
1389- double const* const* const* const* const* const* twoElecTwoCore,
1389+ double const* const* const* const* const* const* twoElecsTwoAtomCores,
13901390 bool isGuess) const{
13911391 int totalNumberAOs = molecule.GetTotalNumberAOs();
13921392 int totalNumberAtoms = molecule.GetNumberAtoms();
@@ -1423,7 +1423,7 @@ void Cndo2::CalcFockMatrix(double** fockMatrix,
14231423 gammaAB,
14241424 orbitalElectronPopulation,
14251425 atomicElectronPopulation,
1426- twoElecTwoCore,
1426+ twoElecsTwoAtomCores,
14271427 isGuess);
14281428 }
14291429 else if(mu < nu){
@@ -1438,7 +1438,7 @@ void Cndo2::CalcFockMatrix(double** fockMatrix,
14381438 gammaAB,
14391439 overlapAOs,
14401440 orbitalElectronPopulation,
1441- twoElecTwoCore,
1441+ twoElecsTwoAtomCores,
14421442 isGuess);
14431443 }
14441444 else{
@@ -1495,7 +1495,7 @@ double Cndo2::GetFockDiagElement(const Atom& atomA,
14951495 double const* const* gammaAB,
14961496 double const* const* orbitalElectronPopulation,
14971497 double const* atomicElectronPopulation,
1498- double const* const* const* const* const* const* twoElecTwoCore,
1498+ double const* const* const* const* const* const* twoElecsTwoAtomCores,
14991499 bool isGuess) const{
15001500 double value;
15011501 int firstAOIndexA = atomA.GetFirstAOIndex();
@@ -1531,7 +1531,7 @@ double Cndo2::GetFockOffDiagElement(const Atom& atomA,
15311531 double const* const* gammaAB,
15321532 double const* const* overlapAOs,
15331533 double const* const* orbitalElectronPopulation,
1534- double const* const* const* const* const* const* twoElecTwoCore,
1534+ double const* const* const* const* const* const* twoElecsTwoAtomCores,
15351535 bool isGuess) const{
15361536 double value;
15371537 double K = this->GetBondingAdjustParameterK(atomA.GetValenceShellType(), atomB.GetValenceShellType());
--- a/src/cndo/Cndo2.h
+++ b/src/cndo/Cndo2.h
@@ -104,7 +104,7 @@ protected:
104104 double** atomicElectronPopulationCIS;
105105 double** atomicUnpairedPopulationCIS;
106106 double** overlapAOs; // overlap integral between AOs
107- double****** twoElecTwoCore;
107+ double****** twoElecsTwoAtomCores;
108108 double*** cartesianMatrix; // cartesian matrix represented by AOs
109109 double*** electronicTransitionDipoleMoments; // Diagnonal terms are electronic dipole moments of each eigenstates (i.e. electronicDipole[0][0][XAxis] is the x-component of the electronic dipole moment of the ground state. electronicDipole[10][10][XAxis] is the x-component of the electronic dipole moment of the 10-th excited state). Off-diagonal terms are transition dipole moments between eigenstates (i.e. electronicDipole[10][0][XAxis] is the x-component of the transition dipole moment from the ground state to 10-th excited state.).
110110 double* coreDipoleMoment; // dipole moment of configuration.
@@ -170,7 +170,7 @@ protected:
170170 double const* const* gammaAB,
171171 double const* const* orbitalElectronPopulation,
172172 double const* atomicElectronPopulation,
173- double const* const* const* const* const* const* twoElecTwoCore,
173+ double const* const* const* const* const* const* twoElecsTwoAtomCores,
174174 bool isGuess) const;
175175 virtual double GetFockOffDiagElement(const MolDS_base_atoms::Atom& atomA,
176176 const MolDS_base_atoms::Atom& atomB,
@@ -182,7 +182,7 @@ protected:
182182 double const* const* gammaAB,
183183 double const* const* overlapAOs,
184184 double const* const* orbitalElectronPopulation,
185- double const* const* const* const* const* const* twoElecTwoCore,
185+ double const* const* const* const* const* const* twoElecsTwoAtomCores,
186186 bool isGuess) const;
187187 void TransposeFockMatrixMatrix(double** transposedFockMatrix) const;
188188 virtual void CalcDiatomicOverlapAOsInDiatomicFrame(double** diatomicOverlapAOs,
@@ -250,8 +250,8 @@ protected:
250250 const MolDS_base::Molecule& molecule,
251251 double const* const* fockMatrix,
252252 double const* const* gammaAB) const;
253- virtual void CalcTwoElecTwoCore(double****** twoElecTwoCore,
254- const MolDS_base::Molecule& molecule) const;
253+ virtual void CalcTwoElecsTwoCores(double****** twoElecsTwoAtomCores,
254+ const MolDS_base::Molecule& molecule) const;
255255 virtual void CalcForce(const std::vector<int>& elecStates);
256256 void CalcRotatingMatrix1stDerivatives(double*** rotMat1stDerivatives,
257257 const MolDS_base_atoms::Atom& atomA,
@@ -463,7 +463,7 @@ private:
463463 double const* const* gammaAB,
464464 double const* const* orbitalElectronPopulation,
465465 double const* atomicElectronPopulation,
466- double const* const* const* const* const* const* twoElecTwoCore,
466+ double const* const* const* const* const* const* twoElecsTwoAtomCores,
467467 bool isGuess) const;
468468 void RotateDiatmicOverlapAOsToSpaceFrame(double** diatomicOverlapAOs,
469469 double const* const* rotatingMatrix,
--- a/src/indo/Indo.cpp
+++ b/src/indo/Indo.cpp
@@ -113,7 +113,7 @@ double Indo::GetFockDiagElement(const Atom& atomA,
113113 double const* const* gammaAB,
114114 double const* const* orbitalElectronPopulation,
115115 double const* atomicElectronPopulation,
116- double const* const* const* const* const* const* twoElecTwoCore,
116+ double const* const* const* const* const* const* twoElecsTwoAtomCores,
117117 bool isGuess) const{
118118 double value;
119119 int firstAOIndexA = atomA.GetFirstAOIndex();
@@ -160,7 +160,7 @@ double Indo::GetFockOffDiagElement(const Atom& atomA,
160160 double const* const* gammaAB,
161161 double const* const* overlapAOs,
162162 double const* const* orbitalElectronPopulation,
163- double const* const* const* const* const* const* twoElecTwoCore,
163+ double const* const* const* const* const* const* twoElecsTwoAtomCores,
164164 bool isGuess) const{
165165 double value;
166166 double K = this->GetBondingAdjustParameterK(atomA.GetValenceShellType(), atomB.GetValenceShellType());
--- a/src/indo/Indo.h
+++ b/src/indo/Indo.h
@@ -37,7 +37,7 @@ protected:
3737 double const* const* gammaAB,
3838 double const* const* orbitalElectronPopulation,
3939 double const* atomicElectronPopulation,
40- double const* const* const* const* const* const* twoElecTwoCore,
40+ double const* const* const* const* const* const* twoElecsTwoAtomCores,
4141 bool isGuess) const;
4242 virtual double GetFockOffDiagElement(const MolDS_base_atoms::Atom& atomA,
4343 const MolDS_base_atoms::Atom& atomB,
@@ -49,7 +49,7 @@ protected:
4949 double const* const* gammaAB,
5050 double const* const* overelap,
5151 double const* const* orbitalElectronPopulation,
52- double const* const* const* const* const* const* twoElecTwoCore,
52+ double const* const* const* const* const* const* twoElecsTwoAtomCores,
5353 bool isGuess) const;
5454 virtual double GetMolecularIntegralElement(int moI, int moJ, int moK, int moL,
5555 const MolDS_base::Molecule& molecule,
--- a/src/mndo/Mndo.cpp
+++ b/src/mndo/Mndo.cpp
@@ -49,6 +49,7 @@
4949 #include"../base/atoms/Natom.h"
5050 #include"../base/atoms/Oatom.h"
5151 #include"../base/atoms/Satom.h"
52+#include"../base/atoms/mm/EnvironmentalPointCharge.h"
5253 #include"../base/Molecule.h"
5354 #include"../base/ElectronicStructure.h"
5455 #include"../cndo/Cndo2.h"
@@ -68,14 +69,14 @@ Mndo::Mndo() : MolDS_zindo::ZindoS(){
6869 this->SetMessages();
6970 this->SetEnableAtomTypes();
7071 // private variables
71- this->twoElecTwoCoreMpiBuff = NULL;
72+ this->twoElecsTwoAtomCoresMpiBuff = NULL;
7273 this->heatsFormation = 0.0;
7374 //this->OutputLog("Mndo created\n");
7475 }
7576
7677 Mndo::~Mndo(){
7778 OrbitalType twoElecLimit = dxy;
78- MallocerFreer::GetInstance()->Free<double>(&this->twoElecTwoCore,
79+ MallocerFreer::GetInstance()->Free<double>(&this->twoElecsTwoAtomCores,
7980 this->molecule->GetNumberAtoms(),
8081 this->molecule->GetNumberAtoms(),
8182 twoElecLimit,
@@ -83,7 +84,7 @@ Mndo::~Mndo(){
8384 twoElecLimit,
8485 twoElecLimit);
8586 int numBuff = (twoElecLimit+1)*twoElecLimit/2;
86- MallocerFreer::GetInstance()->Free<double>(&this->twoElecTwoCoreMpiBuff,
87+ MallocerFreer::GetInstance()->Free<double>(&this->twoElecsTwoAtomCoresMpiBuff,
8788 this->molecule->GetNumberAtoms(),
8889 this->molecule->GetNumberAtoms(),
8990 numBuff,
@@ -98,7 +99,7 @@ Mndo::~Mndo(){
9899 void Mndo::SetMolecule(Molecule* molecule){
99100 ZindoS::SetMolecule(molecule);
100101 OrbitalType twoElecLimit = dxy;
101- MallocerFreer::GetInstance()->Malloc<double>(&this->twoElecTwoCore,
102+ MallocerFreer::GetInstance()->Malloc<double>(&this->twoElecsTwoAtomCores,
102103 molecule->GetNumberAtoms(),
103104 molecule->GetNumberAtoms(),
104105 twoElecLimit,
@@ -106,7 +107,7 @@ void Mndo::SetMolecule(Molecule* molecule){
106107 twoElecLimit,
107108 twoElecLimit);
108109 int numBuff = (twoElecLimit+1)*twoElecLimit/2;
109- MallocerFreer::GetInstance()->Malloc<double>(&this->twoElecTwoCoreMpiBuff,
110+ MallocerFreer::GetInstance()->Malloc<double>(&this->twoElecsTwoAtomCoresMpiBuff,
110111 this->molecule->GetNumberAtoms(),
111112 this->molecule->GetNumberAtoms(),
112113 numBuff,
@@ -145,20 +146,20 @@ void Mndo::SetMessages(){
145146 = "Error in mndo::Mndo::GetNddoRepulsionIntegral1stDerivative: Bad orbital is set.\n";
146147 this->errorMessageGetNddoRepulsionIntegral2ndDerivative
147148 = "Error in mndo::Mndo::GetNddoRepulsionIntegral2ndDerivative: Bad orbital is set.\n";
148- this->errorMessageCalcTwoElecTwoCoreNullMatrix
149- = "Error in mndo::Mndo::CalcTwoElecTwoCore: The two elec two core matrix is NULL.\n";
150- this->errorMessageCalcDiatomicTwoElecTwoCoreSameAtoms
151- = "Error in mndo::Mndo::CalcDiatomicTwoElecTwoCore: Atom A and B is same.\n";
152- this->errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesSameAtoms
153- = "Error in mndo::Mndo::CalcDiatomicTwoElecTwoCore1stDerivatives: Atom A and B is same.\n";
154- this->errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesSameAtoms
155- = "Error in mndo::Mndo::CalcDiatomicTwoElecTwoCore2ndDerivatives: Atom A and B is same.\n";
156- this->errorMessageCalcDiatomicTwoElecTwoCoreNullMatrix
157- = "Error in mndo::Mndo::CalcDiatomicTwoElecTwoCore: The two elec two core diatomic matrix is NULL.\n";
158- this->errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesNullMatrix
159- = "Error in mndo::Mndo::CalcDiatomicTwoElecTwoCore1stDerivatives: The two elec two core diatomic matrix is NULL.\n";
160- this->errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesNullMatrix
161- = "Error in mndo::Mndo::CalcDiatomicTwoElecTwoCore2ndDerivatives: The two elec two core diatomic matrix is NULL.\n";
149+ this->errorMessageCalcTwoElecsTwoCoresNullMatrix
150+ = "Error in mndo::Mndo::CalcTwoElecsTwoCores: The two elec two core matrix is NULL.\n";
151+ this->errorMessageCalcDiatomicTwoElecsTwoCoresSameAtoms
152+ = "Error in mndo::Mndo::CalcDiatomicTwoElecsTwoCores: Atom A and B is same.\n";
153+ this->errorMessageCalcDiatomicTwoElecsTwoCores1stDerivativesSameAtoms
154+ = "Error in mndo::Mndo::CalcDiatomicTwoElecsTwoCores1stDerivatives: Atom A and B is same.\n";
155+ this->errorMessageCalcDiatomicTwoElecsTwoCores2ndDerivativesSameAtoms
156+ = "Error in mndo::Mndo::CalcDiatomicTwoElecsTwoCores2ndDerivatives: Atom A and B is same.\n";
157+ this->errorMessageCalcDiatomicTwoElecsTwoCoresNullMatrix
158+ = "Error in mndo::Mndo::CalcDiatomicTwoElecsTwoCores: The two elec two core diatomic matrix is NULL.\n";
159+ this->errorMessageCalcDiatomicTwoElecsTwoCores1stDerivativesNullMatrix
160+ = "Error in mndo::Mndo::CalcDiatomicTwoElecsTwoCores1stDerivatives: The two elec two core diatomic matrix is NULL.\n";
161+ this->errorMessageCalcDiatomicTwoElecsTwoCores2ndDerivativesNullMatrix
162+ = "Error in mndo::Mndo::CalcDiatomicTwoElecsTwoCores2ndDerivatives: The two elec two core diatomic matrix is NULL.\n";
162163 this->errorMessageGetElectronicEnergyEnergyNotCalculated
163164 = "Error in mndo::Mndo::GetElectronicEnergy: Set electronic state is not calculated by CIS.\n";
164165 this->errorMessageGetElectronicEnergyNULLCISEnergy
@@ -402,7 +403,7 @@ double Mndo::GetFockDiagElement(const Atom& atomA,
402403 double const* const* gammaAB,
403404 double const* const* orbitalElectronPopulation,
404405 double const* atomicElectronPopulation,
405- double const* const* const* const* const* const* twoElecTwoCore,
406+ double const* const* const* const* const* const* twoElecsTwoAtomCores,
406407 bool isGuess) const{
407408 double value=0.0;
408409 int firstAOIndexA = atomA.GetFirstAOIndex();
@@ -432,22 +433,56 @@ double Mndo::GetFockDiagElement(const Atom& atomA,
432433 for(int sigma=0; sigma<valenceSizeB; sigma++){
433434 temp += orbitalElectronPopulation[lambda+firstAOIndexB]
434435 [sigma+firstAOIndexB]
435- *twoElecTwoCore[indexAtomA][B][mu][mu][lambda][sigma];
436+ *twoElecsTwoAtomCores[indexAtomA][B][mu][mu][lambda][sigma];
436437 }
437438 /*
438439 temp += MolDS_wrappers::Blas::GetInstance()->Ddot(valenceSizeB,
439440 &orbitalElectronPopulation[lambda+firstAOIndexB][firstAOIndexB],
440- &twoElecTwoCore[indexAtomA][B][mu][mu][lambda][0]);
441+ &twoElecsTwoAtomCores[indexAtomA][B][mu][mu][lambda][0]);
441442 */
442443 }
443444 temp += this->GetElectronCoreAttraction(indexAtomA,
444445 B,
445446 mu,
446447 mu,
447- twoElecTwoCore);
448+ twoElecsTwoAtomCores);
448449 }
449450 }
450451 value += temp;
452+
453+ /* coulomb repulsion with point charge *
454+ {
455+ Atom* pointCharge = new MolDS_base_atoms_mm::EnvironmentalPointCharge(1000);
456+ pointCharge->SetXyz(0.0,0.0,0.0);
457+ pointCharge->SetPxyz(0.0,0.0,0.0);
458+
459+ double**** diatomicTwoElecsTwoCores = NULL;
460+ double* tmpDiatomicTwoElecsTwoCores = NULL;
461+ double** tmpRotMat = NULL;
462+ double** tmpMatrixBC = NULL;
463+ double* tmpVectorBC = NULL;
464+ MallocerFreer::GetInstance()->Malloc<double>(&diatomicTwoElecsTwoCores, dxy, dxy, dxy, dxy);
465+ MallocerFreer::GetInstance()->Malloc<double>(&tmpDiatomicTwoElecsTwoCores, dxy*dxy*dxy*dxy);
466+ MallocerFreer::GetInstance()->Malloc<double>(&tmpRotMat, OrbitalType_end, OrbitalType_end);
467+ MallocerFreer::GetInstance()->Malloc<double>(&tmpMatrixBC, dxy*dxy, dxy*dxy);
468+ MallocerFreer::GetInstance()->Malloc<double>(&tmpVectorBC, dxy*dxy*dxy*dxy);
469+ this->CalcDiatomicTwoElecsTwoCoresPointCharge(diatomicTwoElecsTwoCores,
470+ tmpDiatomicTwoElecsTwoCores,
471+ tmpRotMat,
472+ tmpMatrixBC,
473+ tmpVectorBC,
474+ atomA,
475+ *pointCharge);
476+ value += diatomicTwoElecsTwoCores[mu][mu][s][s];
477+ MallocerFreer::GetInstance()->Free<double>(&diatomicTwoElecsTwoCores, dxy, dxy, dxy, dxy);
478+ MallocerFreer::GetInstance()->Free<double>(&tmpDiatomicTwoElecsTwoCores, dxy*dxy*dxy*dxy);
479+ MallocerFreer::GetInstance()->Free<double>(&tmpRotMat, OrbitalType_end, OrbitalType_end);
480+ MallocerFreer::GetInstance()->Free<double>(&tmpMatrixBC, dxy*dxy, dxy*dxy);
481+ MallocerFreer::GetInstance()->Free<double>(&tmpVectorBC, dxy*dxy*dxy*dxy);
482+
483+ delete pointCharge;
484+ }
485+ */
451486 }
452487 return value;
453488 }
@@ -462,7 +497,7 @@ double Mndo::GetFockOffDiagElement(const Atom& atomA,
462497 double const* const* gammaAB,
463498 double const* const* overlapAOs,
464499 double const* const* orbitalElectronPopulation,
465- double const* const* const* const* const* const* twoElecTwoCore,
500+ double const* const* const* const* const* const* twoElecsTwoAtomCores,
466501 bool isGuess) const{
467502 double value = 0.0;
468503 int firstAOIndexA = atomA.GetFirstAOIndex();
@@ -495,21 +530,54 @@ double Mndo::GetFockOffDiagElement(const Atom& atomA,
495530 for(int sigma=0; sigma<valenceSizeBB; sigma++){
496531 temp += orbitalElectronPopulation[lambda+firstAOIndexBB]
497532 [sigma+firstAOIndexBB]
498- *twoElecTwoCore[indexAtomA][BB][mu][nu][lambda][sigma];
533+ *twoElecsTwoAtomCores[indexAtomA][BB][mu][nu][lambda][sigma];
499534 }
500535 /*
501536 temp += MolDS_wrappers::Blas::GetInstance()->Ddot(valenceSizeBB,
502537 &orbitalElectronPopulation[lambda+firstAOIndexBB][firstAOIndexBB],
503- &twoElecTwoCore[indexAtomA][BB][mu][nu][lambda][0]);
538+ &twoElecsTwoAtomCores[indexAtomA][BB][mu][nu][lambda][0]);
504539 */
505540 }
506541 temp += this->GetElectronCoreAttraction(indexAtomA,
507542 BB,
508543 mu,
509544 nu,
510- twoElecTwoCore);
545+ twoElecsTwoAtomCores);
511546 }
512547 }
548+ /* coulomb repulsion with point charge *
549+ {
550+ Atom* pointCharge = new MolDS_base_atoms_mm::EnvironmentalPointCharge(1000);
551+ pointCharge->SetXyz(0.0,0.0,0.0);
552+ pointCharge->SetPxyz(0.0,0.0,0.0);
553+
554+ double**** diatomicTwoElecsTwoCores = NULL;
555+ double* tmpDiatomicTwoElecsTwoCores = NULL;
556+ double** tmpRotMat = NULL;
557+ double** tmpMatrixBC = NULL;
558+ double* tmpVectorBC = NULL;
559+ MallocerFreer::GetInstance()->Malloc<double>(&diatomicTwoElecsTwoCores, dxy, dxy, dxy, dxy);
560+ MallocerFreer::GetInstance()->Malloc<double>(&tmpDiatomicTwoElecsTwoCores, dxy*dxy*dxy*dxy);
561+ MallocerFreer::GetInstance()->Malloc<double>(&tmpRotMat, OrbitalType_end, OrbitalType_end);
562+ MallocerFreer::GetInstance()->Malloc<double>(&tmpMatrixBC, dxy*dxy, dxy*dxy);
563+ MallocerFreer::GetInstance()->Malloc<double>(&tmpVectorBC, dxy*dxy*dxy*dxy);
564+ this->CalcDiatomicTwoElecsTwoCoresPointCharge(diatomicTwoElecsTwoCores,
565+ tmpDiatomicTwoElecsTwoCores,
566+ tmpRotMat,
567+ tmpMatrixBC,
568+ tmpVectorBC,
569+ atomA,
570+ *pointCharge);
571+ value += diatomicTwoElecsTwoCores[mu][nu][s][s];
572+ MallocerFreer::GetInstance()->Free<double>(&diatomicTwoElecsTwoCores, dxy, dxy, dxy, dxy);
573+ MallocerFreer::GetInstance()->Free<double>(&tmpDiatomicTwoElecsTwoCores, dxy*dxy*dxy*dxy);
574+ MallocerFreer::GetInstance()->Free<double>(&tmpRotMat, OrbitalType_end, OrbitalType_end);
575+ MallocerFreer::GetInstance()->Free<double>(&tmpMatrixBC, dxy*dxy, dxy*dxy);
576+ MallocerFreer::GetInstance()->Free<double>(&tmpVectorBC, dxy*dxy*dxy*dxy);
577+
578+ delete pointCharge;
579+ }
580+ */
513581 }
514582 else{
515583 temp = bondParameter*overlapAOs[mu+firstAOIndexA][nu+firstAOIndexB];
@@ -518,15 +586,15 @@ double Mndo::GetFockOffDiagElement(const Atom& atomA,
518586 for(int lambda=0; lambda<valenceSizeB; lambda++){
519587 //temp -= 0.5*orbitalElectronPopulation[lambda+firstAOIndexB]
520588 // [sigma+firstAOIndexA]
521- // *twoElecTwoCore[indexAtomA][indexAtomB][mu][sigma][nu][lambda];
589+ // *twoElecsTwoAtomCores[indexAtomA][indexAtomB][mu][sigma][nu][lambda];
522590 temp -= 0.5*orbitalElectronPopulation[sigma+firstAOIndexA]
523591 [lambda+firstAOIndexB]
524- *twoElecTwoCore[indexAtomA][indexAtomB][mu][sigma][nu][lambda];
592+ *twoElecsTwoAtomCores[indexAtomA][indexAtomB][mu][sigma][nu][lambda];
525593 }
526594 /*
527595 temp -= 0.5*MolDS_wrappers::Blas::GetInstance()->Ddot(valenceSizeB,
528596 &orbitalElectronPopulation[sigma+firstAOIndexA][firstAOIndexB],
529- &twoElecTwoCore[indexAtomA][indexAtomB][mu][sigma][nu][0]);
597+ &twoElecsTwoAtomCores[indexAtomA][indexAtomB][mu][sigma][nu][0]);
530598 */
531599 }
532600 }
@@ -600,24 +668,24 @@ double Mndo::GetElectronCoreAttraction(int indexAtomA,
600668 int indexAtomB,
601669 int mu,
602670 int nu,
603- double const* const* const* const* const* const* twoElecTwoCore) const{
671+ double const* const* const* const* const* const* twoElecsTwoAtomCores) const{
604672 const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
605- return -1.0*atomB.GetCoreCharge()*twoElecTwoCore[indexAtomA][indexAtomB][mu][nu][s][s];
673+ return -1.0*atomB.GetCoreCharge()*twoElecsTwoAtomCores[indexAtomA][indexAtomB][mu][nu][s][s];
606674 }
607675
608676 // First derivative of electron in atom A (mu and nu) and core (atom B) attraction.
609677 // This derivative is related to the coordinate of atomA.
610-// Note that diatomicTwoElecTwoCore1stDerivative is dioatomic one.
678+// Note that diatomicTwoElecsTwoCores1stDerivative is dioatomic one.
611679 // see Eq. (16) in [DT_1977-2] with f_2 = 0.
612680 double Mndo::GetElectronCoreAttraction1stDerivative(int indexAtomA,
613681 int indexAtomB,
614682 int mu,
615683 int nu,
616- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivatives,
684+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivatives,
617685 CartesianType axisA) const{
618686 const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
619687 double value = -1.0*atomB.GetCoreCharge()
620- *diatomicTwoElecTwoCore1stDerivatives[mu][nu][s][s][axisA];
688+ *diatomicTwoElecsTwoCores1stDerivatives[mu][nu][s][s][axisA];
621689 return value;
622690 }
623691
@@ -666,12 +734,12 @@ double Mndo::GetMolecularIntegralElement(int moI, int moJ, int moK, int moL,
666734 for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
667735 for(int sigma=lambda; sigma<=lastAOIndexB; sigma++){
668736 OrbitalType orbitalSigma = atomB.GetValence(sigma-firstAOIndexB);
669- gamma = this->twoElecTwoCore[A]
670- [B]
671- [mu-firstAOIndexA]
672- [nu-firstAOIndexA]
673- [lambda-firstAOIndexB]
674- [sigma-firstAOIndexB];
737+ gamma = this->twoElecsTwoAtomCores[A]
738+ [B]
739+ [mu-firstAOIndexA]
740+ [nu-firstAOIndexA]
741+ [lambda-firstAOIndexB]
742+ [sigma-firstAOIndexB];
675743
676744 value += gamma*fockMatrix[moI][mu]
677745 *fockMatrix[moJ][nu]
@@ -821,12 +889,12 @@ void Mndo::CalcCISMatrix(double** matrixCIS) const{
821889 double tmpMuNuLamda16 = tmpMuNu12*fockMatrix[moB][lambda];
822890 for(int sigma=lambda; sigma<=lastAOIndexB; sigma++){
823891 OrbitalType orbitalSigma = atomB.GetValence(sigma-firstAOIndexB);
824- gamma = this->twoElecTwoCore[A]
825- [B]
826- [mu-firstAOIndexA]
827- [nu-firstAOIndexA]
828- [lambda-firstAOIndexB]
829- [sigma-firstAOIndexB];
892+ gamma = this->twoElecsTwoAtomCores[A]
893+ [B]
894+ [mu-firstAOIndexA]
895+ [nu-firstAOIndexA]
896+ [lambda-firstAOIndexB]
897+ [sigma-firstAOIndexB];
830898
831899 value += gamma*tmpMuNuLamda01*fockMatrix[moB][sigma];
832900 value += gamma*tmpMuNuLamda02*fockMatrix[moI][sigma];
@@ -993,14 +1061,14 @@ double Mndo::GetCISCoefficientTwoElecIntegral(int k,
9931061
9941062 void Mndo::MallocTempMatricesEachThreadCalcHessianSCF(double***** diatomicOverlapAOs1stDerivs,
9951063 double****** diatomicOverlapAOs2ndDerivs,
996- double******* diatomicTwoElecTwoCore1stDerivs,
997- double******** diatomicTwoElecTwoCore2ndDerivs,
1064+ double******* diatomicTwoElecsTwoCores1stDerivs,
1065+ double******** diatomicTwoElecsTwoCores2ndDerivs,
9981066 double*** tmpRotMat,
9991067 double*** tmpRotMat1stDeriv,
10001068 double**** tmpRotMat1stDerivs,
10011069 double***** tmpRotMat2ndDerivs,
1002- double***** tmpDiatomicTwoElecTwoCore,
1003- double****** tmpDiatomicTwoElecTwoCore1stDerivs,
1070+ double***** tmpDiatomicTwoElecsTwoCores,
1071+ double****** tmpDiatomicTwoElecsTwoCores1stDerivs,
10041072 double*** tmpDiaOverlapAOsInDiaFrame,
10051073 double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
10061074 double*** tmpDiaOverlapAOs2ndDerivInDiaFrame,
@@ -1021,14 +1089,14 @@ void Mndo::MallocTempMatricesEachThreadCalcHessianSCF(double***** diatomicOve
10211089 OrbitalType_end,
10221090 CartesianType_end,
10231091 CartesianType_end);
1024- MallocerFreer::GetInstance()->Malloc<double>(diatomicTwoElecTwoCore1stDerivs,
1092+ MallocerFreer::GetInstance()->Malloc<double>(diatomicTwoElecsTwoCores1stDerivs,
10251093 this->molecule->GetNumberAtoms(),
10261094 dxy,
10271095 dxy,
10281096 dxy,
10291097 dxy,
10301098 CartesianType_end);
1031- MallocerFreer::GetInstance()->Malloc<double>(diatomicTwoElecTwoCore2ndDerivs,
1099+ MallocerFreer::GetInstance()->Malloc<double>(diatomicTwoElecsTwoCores2ndDerivs,
10321100 this->molecule->GetNumberAtoms(),
10331101 dxy,
10341102 dxy,
@@ -1051,12 +1119,12 @@ void Mndo::MallocTempMatricesEachThreadCalcHessianSCF(double***** diatomicOve
10511119 OrbitalType_end,
10521120 CartesianType_end,
10531121 CartesianType_end);
1054- MallocerFreer::GetInstance()->Malloc<double>(tmpDiatomicTwoElecTwoCore,
1122+ MallocerFreer::GetInstance()->Malloc<double>(tmpDiatomicTwoElecsTwoCores,
10551123 dxy,
10561124 dxy,
10571125 dxy,
10581126 dxy);
1059- MallocerFreer::GetInstance()->Malloc<double>(tmpDiatomicTwoElecTwoCore1stDerivs,
1127+ MallocerFreer::GetInstance()->Malloc<double>(tmpDiatomicTwoElecsTwoCores1stDerivs,
10601128 dxy,
10611129 dxy,
10621130 dxy,
@@ -1093,14 +1161,14 @@ void Mndo::MallocTempMatricesEachThreadCalcHessianSCF(double***** diatomicOve
10931161
10941162 void Mndo::FreeTempMatricesEachThreadCalcHessianSCF(double***** diatomicOverlapAOs1stDerivs,
10951163 double****** diatomicOverlapAOs2ndDerivs,
1096- double******* diatomicTwoElecTwoCore1stDerivs,
1097- double******** diatomicTwoElecTwoCore2ndDerivs,
1164+ double******* diatomicTwoElecsTwoCores1stDerivs,
1165+ double******** diatomicTwoElecsTwoCores2ndDerivs,
10981166 double*** tmpRotMat,
10991167 double*** tmpRotMat1stDeriv,
11001168 double**** tmpRotMat1stDerivs,
11011169 double***** tmpRotMat2ndDerivs,
1102- double***** tmpDiatomicTwoElecTwoCore,
1103- double****** tmpDiatomicTwoElecTwoCore1stDerivs,
1170+ double***** tmpDiatomicTwoElecsTwoCores,
1171+ double****** tmpDiatomicTwoElecsTwoCores1stDerivs,
11041172 double*** tmpDiaOverlapAOsInDiaFrame,
11051173 double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
11061174 double*** tmpDiaOverlapAOs2ndDerivInDiaFrame,
@@ -1121,14 +1189,14 @@ void Mndo::FreeTempMatricesEachThreadCalcHessianSCF(double***** diatomicOverl
11211189 OrbitalType_end,
11221190 CartesianType_end,
11231191 CartesianType_end);
1124- MallocerFreer::GetInstance()->Free<double>(diatomicTwoElecTwoCore1stDerivs,
1192+ MallocerFreer::GetInstance()->Free<double>(diatomicTwoElecsTwoCores1stDerivs,
11251193 this->molecule->GetNumberAtoms(),
11261194 dxy,
11271195 dxy,
11281196 dxy,
11291197 dxy,
11301198 CartesianType_end);
1131- MallocerFreer::GetInstance()->Free<double>(diatomicTwoElecTwoCore2ndDerivs,
1199+ MallocerFreer::GetInstance()->Free<double>(diatomicTwoElecsTwoCores2ndDerivs,
11321200 this->molecule->GetNumberAtoms(),
11331201 dxy,
11341202 dxy,
@@ -1151,12 +1219,12 @@ void Mndo::FreeTempMatricesEachThreadCalcHessianSCF(double***** diatomicOverl
11511219 OrbitalType_end,
11521220 CartesianType_end,
11531221 CartesianType_end);
1154- MallocerFreer::GetInstance()->Free<double>(tmpDiatomicTwoElecTwoCore,
1222+ MallocerFreer::GetInstance()->Free<double>(tmpDiatomicTwoElecsTwoCores,
11551223 dxy,
11561224 dxy,
11571225 dxy,
11581226 dxy);
1159- MallocerFreer::GetInstance()->Free<double>(tmpDiatomicTwoElecTwoCore1stDerivs,
1227+ MallocerFreer::GetInstance()->Free<double>(tmpDiatomicTwoElecsTwoCores1stDerivs,
11601228 dxy,
11611229 dxy,
11621230 dxy,
@@ -1202,18 +1270,18 @@ double Mndo::GetAuxiliaryHessianElement1(int mu,
12021270 CartesianType axisA1,
12031271 CartesianType axisA2,
12041272 double const* const* orbitalElectronPopulation,
1205- double const* const* const* const* const* const* diatomicTwoElecTwoCore2ndDerivs) const{
1273+ double const* const* const* const* const* const* diatomicTwoElecsTwoCores2ndDerivs) const{
12061274 const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
12071275 const Atom& atomC = *this->molecule->GetAtom(indexAtomC);
12081276 int firstAOIndexA = atomA.GetFirstAOIndex();
12091277 double value = orbitalElectronPopulation[mu]
12101278 [nu]
1211- *diatomicTwoElecTwoCore2ndDerivs[mu-firstAOIndexA]
1212- [nu-firstAOIndexA]
1213- [s]
1214- [s]
1215- [axisA1]
1216- [axisA2];
1279+ *diatomicTwoElecsTwoCores2ndDerivs[mu-firstAOIndexA]
1280+ [nu-firstAOIndexA]
1281+ [s]
1282+ [s]
1283+ [axisA1]
1284+ [axisA2];
12171285 return value*atomC.GetCoreCharge();
12181286 }
12191287
@@ -1228,7 +1296,7 @@ double Mndo::GetAuxiliaryHessianElement2(int mu,
12281296 CartesianType axisA,
12291297 CartesianType axisB,
12301298 double const* const* const* const* orbitalElectronPopulation1stDerivs,
1231- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs) const{
1299+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
12321300 const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
12331301 const Atom& atomC = *this->molecule->GetAtom(indexAtomC);
12341302 int firstAOIndexA = atomA.GetFirstAOIndex();
@@ -1236,11 +1304,11 @@ double Mndo::GetAuxiliaryHessianElement2(int mu,
12361304 [nu]
12371305 [indexAtomB]
12381306 [axisB]
1239- *diatomicTwoElecTwoCore1stDerivs[mu-firstAOIndexA]
1240- [nu-firstAOIndexA]
1241- [s]
1242- [s]
1243- [axisA];
1307+ *diatomicTwoElecsTwoCores1stDerivs[mu-firstAOIndexA]
1308+ [nu-firstAOIndexA]
1309+ [s]
1310+ [s]
1311+ [axisA];
12441312 return value*atomC.GetCoreCharge();
12451313 }
12461314
@@ -1254,18 +1322,18 @@ double Mndo::GetAuxiliaryHessianElement3(int lambda,
12541322 CartesianType axisA1,
12551323 CartesianType axisA2,
12561324 double const* const* orbitalElectronPopulation,
1257- double const* const* const* const* const* const* diatomicTwoElecTwoCore2ndDerivs) const{
1325+ double const* const* const* const* const* const* diatomicTwoElecsTwoCores2ndDerivs) const{
12581326 const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
12591327 const Atom& atomC = *this->molecule->GetAtom(indexAtomC);
12601328 int firstAOIndexC = atomC.GetFirstAOIndex();
12611329 double value = orbitalElectronPopulation[lambda]
12621330 [sigma]
1263- *diatomicTwoElecTwoCore2ndDerivs[s]
1264- [s]
1265- [lambda-firstAOIndexC]
1266- [sigma-firstAOIndexC]
1267- [axisA1]
1268- [axisA2];
1331+ *diatomicTwoElecsTwoCores2ndDerivs[s]
1332+ [s]
1333+ [lambda-firstAOIndexC]
1334+ [sigma-firstAOIndexC]
1335+ [axisA1]
1336+ [axisA2];
12691337 return value*atomA.GetCoreCharge();
12701338 }
12711339
@@ -1280,7 +1348,7 @@ double Mndo::GetAuxiliaryHessianElement4(int lambda,
12801348 CartesianType axisA,
12811349 CartesianType axisB,
12821350 double const* const* const* const* orbitalElectronPopulation1stDerivs,
1283- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs) const{
1351+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
12841352 const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
12851353 const Atom& atomC = *this->molecule->GetAtom(indexAtomC);
12861354 int firstAOIndexC = atomC.GetFirstAOIndex();
@@ -1288,11 +1356,11 @@ double Mndo::GetAuxiliaryHessianElement4(int lambda,
12881356 [sigma]
12891357 [indexAtomB]
12901358 [axisB]
1291- *diatomicTwoElecTwoCore1stDerivs[s]
1292- [s]
1293- [lambda-firstAOIndexC]
1294- [sigma-firstAOIndexC]
1295- [axisA];
1359+ *diatomicTwoElecsTwoCores1stDerivs[s]
1360+ [s]
1361+ [lambda-firstAOIndexC]
1362+ [sigma-firstAOIndexC]
1363+ [axisA];
12961364 return value*atomA.GetCoreCharge();
12971365 }
12981366
@@ -1365,7 +1433,7 @@ double Mndo::GetAuxiliaryHessianElement7(int mu,
13651433 CartesianType axisA1,
13661434 CartesianType axisA2,
13671435 double const* const* orbitalElectronPopulation,
1368- double const* const* const* const* const* const* diatomicTwoElecTwoCore2ndDerivs) const{
1436+ double const* const* const* const* const* const* diatomicTwoElecsTwoCores2ndDerivs) const{
13691437 const Atom& atomA = *molecule->GetAtom(indexAtomA);
13701438 const Atom& atomC = *molecule->GetAtom(indexAtomC);
13711439 int firstAOIndexA = atomA.GetFirstAOIndex();
@@ -1373,12 +1441,12 @@ double Mndo::GetAuxiliaryHessianElement7(int mu,
13731441 double temp1 = orbitalElectronPopulation[mu][nu]*orbitalElectronPopulation[lambda][sigma];
13741442 double temp2 = orbitalElectronPopulation[mu][lambda]*orbitalElectronPopulation[nu][sigma];
13751443 double value = (temp1 - 0.5*temp2)
1376- *diatomicTwoElecTwoCore2ndDerivs[mu-firstAOIndexA]
1377- [nu-firstAOIndexA]
1378- [lambda-firstAOIndexC]
1379- [sigma-firstAOIndexC]
1380- [axisA1]
1381- [axisA2];
1444+ *diatomicTwoElecsTwoCores2ndDerivs[mu-firstAOIndexA]
1445+ [nu-firstAOIndexA]
1446+ [lambda-firstAOIndexC]
1447+ [sigma-firstAOIndexC]
1448+ [axisA1]
1449+ [axisA2];
13821450 return value;
13831451 }
13841452
@@ -1396,7 +1464,7 @@ double Mndo::GetAuxiliaryHessianElement8(int mu,
13961464 CartesianType axisB,
13971465 double const* const* orbitalElectronPopulation,
13981466 double const* const* const* const* orbitalElectronPopulation1stDerivs,
1399- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs) const{
1467+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
14001468 const Atom& atomA = *molecule->GetAtom(indexAtomA);
14011469 const Atom& atomC = *molecule->GetAtom(indexAtomC);
14021470 int firstAOIndexA = atomA.GetFirstAOIndex();
@@ -1410,11 +1478,11 @@ double Mndo::GetAuxiliaryHessianElement8(int mu,
14101478 double temp4 = orbitalElectronPopulation [mu][lambda]
14111479 *orbitalElectronPopulation1stDerivs[nu][sigma] [indexAtomB][axisB];
14121480 double value = ((temp1 + temp2) - 0.5*(temp3 + temp4))
1413- *diatomicTwoElecTwoCore1stDerivs[mu-firstAOIndexA]
1414- [nu-firstAOIndexA]
1415- [lambda-firstAOIndexC]
1416- [sigma-firstAOIndexC]
1417- [axisA];
1481+ *diatomicTwoElecsTwoCores1stDerivs[mu-firstAOIndexA]
1482+ [nu-firstAOIndexA]
1483+ [lambda-firstAOIndexC]
1484+ [sigma-firstAOIndexC]
1485+ [axisA];
14181486 return value;
14191487 }
14201488
@@ -1428,8 +1496,8 @@ double Mndo::GetHessianElementSameAtomsSCF(int indexAtomA,
14281496 double const* const* const* const* orbitalElectronPopulation1stDerivs,
14291497 double const* const* const* const* diatomicOverlapAOs1stDerivs,
14301498 double const* const* const* const* const* diatomicOverlapAOs2ndDerivs,
1431- double const* const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs,
1432- double const* const* const* const* const* const* const* diatomicTwoElecTwoCore2ndDerivs) const{
1499+ double const* const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs,
1500+ double const* const* const* const* const* const* const* diatomicTwoElecsTwoCores2ndDerivs) const{
14331501 double value=0.0;
14341502 int indexAtomB = indexAtomA;
14351503 const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
@@ -1451,7 +1519,7 @@ double Mndo::GetHessianElementSameAtomsSCF(int indexAtomA,
14511519 static_cast<CartesianType>(axisA1),
14521520 static_cast<CartesianType>(axisA2),
14531521 orbitalElectronPopulation,
1454- diatomicTwoElecTwoCore2ndDerivs[indexAtomC]);
1522+ diatomicTwoElecsTwoCores2ndDerivs[indexAtomC]);
14551523 value -= this->GetAuxiliaryHessianElement2(mu,
14561524 nu,
14571525 indexAtomA,
@@ -1460,7 +1528,7 @@ double Mndo::GetHessianElementSameAtomsSCF(int indexAtomA,
14601528 static_cast<CartesianType>(axisA1),
14611529 static_cast<CartesianType>(axisA2),
14621530 orbitalElectronPopulation1stDerivs,
1463- diatomicTwoElecTwoCore1stDerivs[indexAtomC]);
1531+ diatomicTwoElecsTwoCores1stDerivs[indexAtomC]);
14641532 }
14651533 }
14661534 for(int lambda=firstAOIndexC; lambda<firstAOIndexC+numberAOsC; lambda++){
@@ -1472,7 +1540,7 @@ double Mndo::GetHessianElementSameAtomsSCF(int indexAtomA,
14721540 static_cast<CartesianType>(axisA1),
14731541 static_cast<CartesianType>(axisA2),
14741542 orbitalElectronPopulation,
1475- diatomicTwoElecTwoCore2ndDerivs[indexAtomC]);
1543+ diatomicTwoElecsTwoCores2ndDerivs[indexAtomC]);
14761544 value -= this->GetAuxiliaryHessianElement4(lambda,
14771545 sigma,
14781546 indexAtomA,
@@ -1481,7 +1549,7 @@ double Mndo::GetHessianElementSameAtomsSCF(int indexAtomA,
14811549 static_cast<CartesianType>(axisA1),
14821550 static_cast<CartesianType>(axisA2),
14831551 orbitalElectronPopulation1stDerivs,
1484- diatomicTwoElecTwoCore1stDerivs[indexAtomC]);
1552+ diatomicTwoElecsTwoCores1stDerivs[indexAtomC]);
14851553 }
14861554 }
14871555 for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
@@ -1518,7 +1586,7 @@ double Mndo::GetHessianElementSameAtomsSCF(int indexAtomA,
15181586 static_cast<CartesianType>(axisA1),
15191587 static_cast<CartesianType>(axisA2),
15201588 orbitalElectronPopulation,
1521- diatomicTwoElecTwoCore2ndDerivs[indexAtomC]);
1589+ diatomicTwoElecsTwoCores2ndDerivs[indexAtomC]);
15221590 value += this->GetAuxiliaryHessianElement8(mu,
15231591 nu,
15241592 lambda,
@@ -1530,7 +1598,7 @@ double Mndo::GetHessianElementSameAtomsSCF(int indexAtomA,
15301598 static_cast<CartesianType>(axisA2),
15311599 orbitalElectronPopulation,
15321600 orbitalElectronPopulation1stDerivs,
1533- diatomicTwoElecTwoCore1stDerivs[indexAtomC]);
1601+ diatomicTwoElecsTwoCores1stDerivs[indexAtomC]);
15341602 }
15351603 }
15361604 }
@@ -1565,8 +1633,8 @@ double Mndo::GetHessianElementDifferentAtomsSCF(int indexAtomA,
15651633 double const* const* const* const* orbitalElectronPopulation1stDerivs,
15661634 double const* const* const* const* diatomicOverlapAOs1stDerivs,
15671635 double const* const* const* const* const* diatomicOverlapAOs2ndDerivs,
1568- double const* const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs,
1569- double const* const* const* const* const* const* const* diatomicTwoElecTwoCore2ndDerivs) const{
1636+ double const* const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs,
1637+ double const* const* const* const* const* const* const* diatomicTwoElecsTwoCores2ndDerivs) const{
15701638 double value=0.0;
15711639 const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
15721640 const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
@@ -1585,7 +1653,7 @@ double Mndo::GetHessianElementDifferentAtomsSCF(int indexAtomA,
15851653 static_cast<CartesianType>(axisA),
15861654 static_cast<CartesianType>(axisB),
15871655 orbitalElectronPopulation,
1588- diatomicTwoElecTwoCore2ndDerivs[indexAtomB]);
1656+ diatomicTwoElecsTwoCores2ndDerivs[indexAtomB]);
15891657 }
15901658 }
15911659 for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
@@ -1597,7 +1665,7 @@ double Mndo::GetHessianElementDifferentAtomsSCF(int indexAtomA,
15971665 static_cast<CartesianType>(axisA),
15981666 static_cast<CartesianType>(axisB),
15991667 orbitalElectronPopulation,
1600- diatomicTwoElecTwoCore2ndDerivs[indexAtomB]);
1668+ diatomicTwoElecsTwoCores2ndDerivs[indexAtomB]);
16011669 }
16021670 }
16031671 for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
@@ -1625,7 +1693,7 @@ double Mndo::GetHessianElementDifferentAtomsSCF(int indexAtomA,
16251693 static_cast<CartesianType>(axisA),
16261694 static_cast<CartesianType>(axisB),
16271695 orbitalElectronPopulation,
1628- diatomicTwoElecTwoCore2ndDerivs[indexAtomB]);
1696+ diatomicTwoElecsTwoCores2ndDerivs[indexAtomB]);
16291697 }
16301698 }
16311699 }
@@ -1648,7 +1716,7 @@ double Mndo::GetHessianElementDifferentAtomsSCF(int indexAtomA,
16481716 static_cast<CartesianType>(axisA),
16491717 static_cast<CartesianType>(axisB),
16501718 orbitalElectronPopulation1stDerivs,
1651- diatomicTwoElecTwoCore1stDerivs[indexAtomC]);
1719+ diatomicTwoElecsTwoCores1stDerivs[indexAtomC]);
16521720 }
16531721 }
16541722 for(int lambda=firstAOIndexC; lambda<firstAOIndexC+numberAOsC; lambda++){
@@ -1661,7 +1729,7 @@ double Mndo::GetHessianElementDifferentAtomsSCF(int indexAtomA,
16611729 static_cast<CartesianType>(axisA),
16621730 static_cast<CartesianType>(axisB),
16631731 orbitalElectronPopulation1stDerivs,
1664- diatomicTwoElecTwoCore1stDerivs[indexAtomC]);
1732+ diatomicTwoElecsTwoCores1stDerivs[indexAtomC]);
16651733 }
16661734 }
16671735 for(int mu=firstAOIndexA; mu<=lastAOIndexA; mu++){
@@ -1692,7 +1760,7 @@ double Mndo::GetHessianElementDifferentAtomsSCF(int indexAtomA,
16921760 static_cast<CartesianType>(axisB),
16931761 orbitalElectronPopulation,
16941762 orbitalElectronPopulation1stDerivs,
1695- diatomicTwoElecTwoCore1stDerivs[indexAtomC]);
1763+ diatomicTwoElecsTwoCores1stDerivs[indexAtomC]);
16961764 }
16971765 }
16981766 }
@@ -1731,37 +1799,37 @@ void Mndo::CalcHessianSCF(double** hessianSCF, bool isMassWeighted) const{
17311799 stringstream ompErrors;
17321800 #pragma omp parallel
17331801 {
1734- double**** diatomicOverlapAOs1stDerivs = NULL;
1735- double***** diatomicOverlapAOs2ndDerivs = NULL;
1736- double****** diatomicTwoElecTwoCore1stDerivs = NULL;
1737- double******* diatomicTwoElecTwoCore2ndDerivs = NULL;
1738- double** tmpRotMat = NULL;
1739- double*** tmpRotMat1stDerivs = NULL;
1740- double**** tmpRotMat2ndDerivs = NULL;
1741- double**** tmpDiatomicTwoElecTwoCore = NULL;
1742- double***** tmpDiatomicTwoElecTwoCore1stDerivs = NULL;
1743- double** tmpDiaOverlapAOsInDiaFrame = NULL; // diatomic overlapAOs in diatomic frame
1744- double** tmpDiaOverlapAOs1stDerivInDiaFrame = NULL; // first derivative of the diaOverlapAOs. This derivative is related to the distance between two atoms.
1745- double** tmpDiaOverlapAOs2ndDerivInDiaFrame = NULL; // second derivative of the diaOverlapAOs. This derivative is related to the distance between two atoms.
1746- double*** tmpDiaOverlapAOs1stDerivs = NULL; // first derivatives of the diaOverlapAOs. This derivatives are related to the all Cartesian coordinates.
1747- double**** tmpDiaOverlapAOs2ndDerivs = NULL; //sedond derivatives of the diaOverlapAOs. This derivatives are related to the all Cartesian coordinates.
1748- double** tmpRotMat1stDeriv = NULL;
1749- double** tmpRotatedDiatomicOverlap = NULL;
1750- double* tmpRotatedDiatomicOverlapVec = NULL; // used in dgemmm
1751- double** tmpMatrixBC = NULL; // used in dgemmm
1752- double* tmpVectorBC = NULL; // used in dgemmm
1802+ double**** diatomicOverlapAOs1stDerivs = NULL;
1803+ double***** diatomicOverlapAOs2ndDerivs = NULL;
1804+ double****** diatomicTwoElecsTwoCores1stDerivs = NULL;
1805+ double******* diatomicTwoElecsTwoCores2ndDerivs = NULL;
1806+ double** tmpRotMat = NULL;
1807+ double*** tmpRotMat1stDerivs = NULL;
1808+ double**** tmpRotMat2ndDerivs = NULL;
1809+ double**** tmpDiatomicTwoElecsTwoCores = NULL;
1810+ double***** tmpDiatomicTwoElecsTwoCores1stDerivs = NULL;
1811+ double** tmpDiaOverlapAOsInDiaFrame = NULL; // diatomic overlapAOs in diatomic frame
1812+ double** tmpDiaOverlapAOs1stDerivInDiaFrame = NULL; // first derivative of the diaOverlapAOs. This derivative is related to the distance between two atoms.
1813+ double** tmpDiaOverlapAOs2ndDerivInDiaFrame = NULL; // second derivative of the diaOverlapAOs. This derivative is related to the distance between two atoms.
1814+ double*** tmpDiaOverlapAOs1stDerivs = NULL; // first derivatives of the diaOverlapAOs. This derivatives are related to the all Cartesian coordinates.
1815+ double**** tmpDiaOverlapAOs2ndDerivs = NULL; //sedond derivatives of the diaOverlapAOs. This derivatives are related to the all Cartesian coordinates.
1816+ double** tmpRotMat1stDeriv = NULL;
1817+ double** tmpRotatedDiatomicOverlap = NULL;
1818+ double* tmpRotatedDiatomicOverlapVec = NULL; // used in dgemmm
1819+ double** tmpMatrixBC = NULL; // used in dgemmm
1820+ double* tmpVectorBC = NULL; // used in dgemmm
17531821
17541822 try{
17551823 this->MallocTempMatricesEachThreadCalcHessianSCF(&diatomicOverlapAOs1stDerivs,
17561824 &diatomicOverlapAOs2ndDerivs,
1757- &diatomicTwoElecTwoCore1stDerivs,
1758- &diatomicTwoElecTwoCore2ndDerivs,
1825+ &diatomicTwoElecsTwoCores1stDerivs,
1826+ &diatomicTwoElecsTwoCores2ndDerivs,
17591827 &tmpRotMat,
17601828 &tmpRotMat1stDeriv,
17611829 &tmpRotMat1stDerivs,
17621830 &tmpRotMat2ndDerivs,
1763- &tmpDiatomicTwoElecTwoCore,
1764- &tmpDiatomicTwoElecTwoCore1stDerivs,
1831+ &tmpDiatomicTwoElecsTwoCores,
1832+ &tmpDiatomicTwoElecsTwoCores1stDerivs,
17651833 &tmpDiaOverlapAOsInDiaFrame,
17661834 &tmpDiaOverlapAOs1stDerivInDiaFrame,
17671835 &tmpDiaOverlapAOs2ndDerivInDiaFrame,
@@ -1804,20 +1872,20 @@ void Mndo::CalcHessianSCF(double** hessianSCF, bool isMassWeighted) const{
18041872 tmpRotMat2ndDerivs,
18051873 indexAtomA,
18061874 indexAtomB);
1807- this->CalcDiatomicTwoElecTwoCore1stDerivatives(diatomicTwoElecTwoCore1stDerivs[indexAtomB],
1808- tmpRotMat,
1809- tmpRotMat1stDerivs,
1810- tmpDiatomicTwoElecTwoCore,
1811- indexAtomA,
1812- indexAtomB);
1813- this->CalcDiatomicTwoElecTwoCore2ndDerivatives(diatomicTwoElecTwoCore2ndDerivs[indexAtomB],
1814- tmpRotMat,
1815- tmpRotMat1stDerivs,
1816- tmpRotMat2ndDerivs,
1817- tmpDiatomicTwoElecTwoCore,
1818- tmpDiatomicTwoElecTwoCore1stDerivs,
1819- indexAtomA,
1820- indexAtomB);
1875+ this->CalcDiatomicTwoElecsTwoCores1stDerivatives(diatomicTwoElecsTwoCores1stDerivs[indexAtomB],
1876+ tmpRotMat,
1877+ tmpRotMat1stDerivs,
1878+ tmpDiatomicTwoElecsTwoCores,
1879+ indexAtomA,
1880+ indexAtomB);
1881+ this->CalcDiatomicTwoElecsTwoCores2ndDerivatives(diatomicTwoElecsTwoCores2ndDerivs[indexAtomB],
1882+ tmpRotMat,
1883+ tmpRotMat1stDerivs,
1884+ tmpRotMat2ndDerivs,
1885+ tmpDiatomicTwoElecsTwoCores,
1886+ tmpDiatomicTwoElecsTwoCores1stDerivs,
1887+ indexAtomA,
1888+ indexAtomB);
18211889 }
18221890 }
18231891
@@ -1837,8 +1905,8 @@ void Mndo::CalcHessianSCF(double** hessianSCF, bool isMassWeighted) const{
18371905 orbitalElectronPopulation1stDerivs,
18381906 diatomicOverlapAOs1stDerivs,
18391907 diatomicOverlapAOs2ndDerivs,
1840- diatomicTwoElecTwoCore1stDerivs,
1841- diatomicTwoElecTwoCore2ndDerivs);
1908+ diatomicTwoElecsTwoCores1stDerivs,
1909+ diatomicTwoElecsTwoCores2ndDerivs);
18421910 if(isMassWeighted){
18431911 hessianSCF[k][l] /= sqrt(atomA.GetCoreMass()*atomB.GetCoreMass());
18441912 }
@@ -1855,8 +1923,8 @@ void Mndo::CalcHessianSCF(double** hessianSCF, bool isMassWeighted) const{
18551923 orbitalElectronPopulation1stDerivs,
18561924 diatomicOverlapAOs1stDerivs,
18571925 diatomicOverlapAOs2ndDerivs,
1858- diatomicTwoElecTwoCore1stDerivs,
1859- diatomicTwoElecTwoCore2ndDerivs);
1926+ diatomicTwoElecsTwoCores1stDerivs,
1927+ diatomicTwoElecsTwoCores2ndDerivs);
18601928 if(isMassWeighted){
18611929 hessianSCF[k][l] /= atomA.GetCoreMass();
18621930 }
@@ -1873,14 +1941,14 @@ void Mndo::CalcHessianSCF(double** hessianSCF, bool isMassWeighted) const{
18731941 }
18741942 this->FreeTempMatricesEachThreadCalcHessianSCF(&diatomicOverlapAOs1stDerivs,
18751943 &diatomicOverlapAOs2ndDerivs,
1876- &diatomicTwoElecTwoCore1stDerivs,
1877- &diatomicTwoElecTwoCore2ndDerivs,
1944+ &diatomicTwoElecsTwoCores1stDerivs,
1945+ &diatomicTwoElecsTwoCores2ndDerivs,
18781946 &tmpRotMat,
18791947 &tmpRotMat1stDeriv,
18801948 &tmpRotMat1stDerivs,
18811949 &tmpRotMat2ndDerivs,
1882- &tmpDiatomicTwoElecTwoCore,
1883- &tmpDiatomicTwoElecTwoCore1stDerivs,
1950+ &tmpDiatomicTwoElecsTwoCores,
1951+ &tmpDiatomicTwoElecsTwoCores1stDerivs,
18841952 &tmpDiaOverlapAOsInDiaFrame,
18851953 &tmpDiaOverlapAOs1stDerivInDiaFrame,
18861954 &tmpDiaOverlapAOs2ndDerivInDiaFrame,
@@ -1988,12 +2056,12 @@ void Mndo::CalcOrbitalElectronPopulation1stDerivatives(double**** orbitalElectro
19882056 for(int indexAtomA=0; indexAtomA<this->molecule->GetNumberAtoms(); indexAtomA++){
19892057 for(int axis=XAxis; axis<CartesianType_end; axis++){
19902058 double temp=0.0;
1991- printf("hoge-cphf: atom:%d axis:%s start\n ",indexAtomA,CartesianTypeStr(axis));
2059+ printf("cphf: atom:%d axis:%s start\n ",indexAtomA,CartesianTypeStr(axis));
19922060 for(int mu=0; mu<totalNumberAOs; mu++){
19932061 temp += orbitalElectronPopulation1stDerivs[mu][mu][indexAtomA][axis];
19942062 printf("%e\n",orbitalElectronPopulation1stDerivs[mu][mu][indexAtomA][axis]);
19952063 }
1996- printf("hoge-cphf: atom:%d axis:%s %e\n\n",indexAtomA,CartesianTypeStr(axis),temp);
2064+ printf("cphf: atom:%d axis:%s %e\n\n",indexAtomA,CartesianTypeStr(axis),temp);
19972065 }
19982066 }
19992067 */
@@ -2070,11 +2138,11 @@ void Mndo::CalcStaticFirstOrderFock(double* staticFirstOrderFock,
20702138 CartesianType axisA) const{
20712139 MallocerFreer::GetInstance()->Initialize<double>(staticFirstOrderFock,
20722140 nonRedundantQIndeces.size()+redundantQIndeces.size());
2073- double***** diatomicTwoElecTwoCore1stDerivs = NULL;
2074- double*** diatomicOverlapAOs1stDerivs = NULL;
2075- double** tmpRotMat = NULL;
2076- double*** tmpRotMat1stDerivs = NULL;
2077- double**** tmpDiatomicTwoElecTwoCore = NULL;
2141+ double***** diatomicTwoElecsTwoCores1stDerivs = NULL;
2142+ double*** diatomicOverlapAOs1stDerivs = NULL;
2143+ double** tmpRotMat = NULL;
2144+ double*** tmpRotMat1stDerivs = NULL;
2145+ double**** tmpDiatomicTwoElecsTwoCores = NULL;
20782146
20792147 double** tmpDiaOverlapAOsInDiaFrame = NULL; // diatomic overlapAOs in diatomic frame
20802148 double** tmpDiaOverlapAOs1stDerivInDiaFrame = NULL; // first derivative of the diaOverlapAOs. This derivative is related to the distance between two atoms.
@@ -2084,11 +2152,11 @@ void Mndo::CalcStaticFirstOrderFock(double* staticFirstOrderFock,
20842152 double** tmpMatrixBC = NULL;
20852153 double* tmpVectorBC = NULL;
20862154 try{
2087- this->MallocTempMatricesStaticFirstOrderFock(&diatomicTwoElecTwoCore1stDerivs,
2155+ this->MallocTempMatricesStaticFirstOrderFock(&diatomicTwoElecsTwoCores1stDerivs,
20882156 &diatomicOverlapAOs1stDerivs,
20892157 &tmpRotMat,
20902158 &tmpRotMat1stDerivs,
2091- &tmpDiatomicTwoElecTwoCore);
2159+ &tmpDiatomicTwoElecsTwoCores);
20922160 MallocerFreer::GetInstance()->Malloc<double>(&tmpDiaOverlapAOsInDiaFrame, OrbitalType_end, OrbitalType_end);
20932161 MallocerFreer::GetInstance()->Malloc<double>(&tmpDiaOverlapAOs1stDerivInDiaFrame, OrbitalType_end, OrbitalType_end);
20942162 MallocerFreer::GetInstance()->Malloc<double>(&tmpRotMat1stDeriv, OrbitalType_end, OrbitalType_end);
@@ -2108,11 +2176,11 @@ void Mndo::CalcStaticFirstOrderFock(double* staticFirstOrderFock,
21082176 int coreChargeB = atomB.GetCoreCharge();
21092177
21102178 // calc. first derivative of two elec two core interaction
2111- this->CalcDiatomicTwoElecTwoCore1stDerivatives(diatomicTwoElecTwoCore1stDerivs,
2112- tmpRotMat,
2113- tmpRotMat1stDerivs,
2114- tmpDiatomicTwoElecTwoCore,
2115- indexAtomA, indexAtomB);
2179+ this->CalcDiatomicTwoElecsTwoCores1stDerivatives(diatomicTwoElecsTwoCores1stDerivs,
2180+ tmpRotMat,
2181+ tmpRotMat1stDerivs,
2182+ tmpDiatomicTwoElecsTwoCores,
2183+ indexAtomA, indexAtomB);
21162184 // calc. first derivative of overlapAOs.
21172185 this->CalcDiatomicOverlapAOs1stDerivatives(diatomicOverlapAOs1stDerivs,
21182186 tmpDiaOverlapAOsInDiaFrame,
@@ -2157,22 +2225,22 @@ void Mndo::CalcStaticFirstOrderFock(double* staticFirstOrderFock,
21572225 *this->fockMatrix[moI][lambda]
21582226 *this->fockMatrix[moJ][mu]
21592227 *this->orbitalElectronPopulation[nu][sigma];
2160- staticFirstOrderFock[i] += temp1*diatomicTwoElecTwoCore1stDerivs[mu-firstAOIndexA]
2161- [nu-firstAOIndexA]
2162- [lambda-firstAOIndexB]
2163- [sigma-firstAOIndexB]
2164- [axisA];
2228+ staticFirstOrderFock[i] += temp1*diatomicTwoElecsTwoCores1stDerivs[mu-firstAOIndexA]
2229+ [nu-firstAOIndexA]
2230+ [lambda-firstAOIndexB]
2231+ [sigma-firstAOIndexB]
2232+ [axisA];
21652233 } //sigma-loop
21662234 } // lambda-loop
21672235
21682236 double temp2 = this->fockMatrix[moI][mu]
21692237 *this->fockMatrix[moJ][nu]
21702238 *coreChargeB
2171- *diatomicTwoElecTwoCore1stDerivs[mu-firstAOIndexA]
2172- [nu-firstAOIndexA]
2173- [s]
2174- [s]
2175- [axisA];
2239+ *diatomicTwoElecsTwoCores1stDerivs[mu-firstAOIndexA]
2240+ [nu-firstAOIndexA]
2241+ [s]
2242+ [s]
2243+ [axisA];
21762244 staticFirstOrderFock[i] -= temp2;
21772245
21782246 } // nu-loop
@@ -2184,11 +2252,11 @@ void Mndo::CalcStaticFirstOrderFock(double* staticFirstOrderFock,
21842252 double temp3 = this->fockMatrix[moI][lambda]
21852253 *this->fockMatrix[moJ][sigma]
21862254 *coreChargeA
2187- *diatomicTwoElecTwoCore1stDerivs[s]
2188- [s]
2189- [lambda-firstAOIndexB]
2190- [sigma-firstAOIndexB]
2191- [axisA];
2255+ *diatomicTwoElecsTwoCores1stDerivs[s]
2256+ [s]
2257+ [lambda-firstAOIndexB]
2258+ [sigma-firstAOIndexB]
2259+ [axisA];
21922260 staticFirstOrderFock[i] -= temp3;
21932261
21942262 } //sigma-loop
@@ -2216,11 +2284,11 @@ void Mndo::CalcStaticFirstOrderFock(double* staticFirstOrderFock,
22162284 }
22172285 }
22182286 catch(MolDSException ex){
2219- this->FreeTempMatricesStaticFirstOrderFock(&diatomicTwoElecTwoCore1stDerivs,
2287+ this->FreeTempMatricesStaticFirstOrderFock(&diatomicTwoElecsTwoCores1stDerivs,
22202288 &diatomicOverlapAOs1stDerivs,
22212289 &tmpRotMat,
22222290 &tmpRotMat1stDerivs,
2223- &tmpDiatomicTwoElecTwoCore);
2291+ &tmpDiatomicTwoElecsTwoCores);
22242292 MallocerFreer::GetInstance()->Free<double>(&tmpDiaOverlapAOsInDiaFrame, OrbitalType_end, OrbitalType_end);
22252293 MallocerFreer::GetInstance()->Free<double>(&tmpDiaOverlapAOs1stDerivInDiaFrame, OrbitalType_end, OrbitalType_end);
22262294 //MallocerFreer::GetInstance()->Free<double>(&tmpRotMat, OrbitalType_end, OrbitalType_end);
@@ -2232,11 +2300,11 @@ void Mndo::CalcStaticFirstOrderFock(double* staticFirstOrderFock,
22322300 MallocerFreer::GetInstance()->Free<double>(&tmpVectorBC, OrbitalType_end*OrbitalType_end);
22332301 throw ex;
22342302 }
2235- this->FreeTempMatricesStaticFirstOrderFock(&diatomicTwoElecTwoCore1stDerivs,
2303+ this->FreeTempMatricesStaticFirstOrderFock(&diatomicTwoElecsTwoCores1stDerivs,
22362304 &diatomicOverlapAOs1stDerivs,
22372305 &tmpRotMat,
22382306 &tmpRotMat1stDerivs,
2239- &tmpDiatomicTwoElecTwoCore);
2307+ &tmpDiatomicTwoElecsTwoCores);
22402308 MallocerFreer::GetInstance()->Free<double>(&tmpDiaOverlapAOsInDiaFrame, OrbitalType_end, OrbitalType_end);
22412309 MallocerFreer::GetInstance()->Free<double>(&tmpDiaOverlapAOs1stDerivInDiaFrame, OrbitalType_end, OrbitalType_end);
22422310 //MallocerFreer::GetInstance()->Free<double>(&tmpRotMat, OrbitalType_end, OrbitalType_end);
@@ -2255,12 +2323,12 @@ void Mndo::CalcStaticFirstOrderFock(double* staticFirstOrderFock,
22552323 */
22562324 }
22572325
2258-void Mndo::MallocTempMatricesStaticFirstOrderFock(double****** diatomicTwoElecTwoCore1stDeriv,
2326+void Mndo::MallocTempMatricesStaticFirstOrderFock(double****** diatomicTwoElecsTwoCores1stDeriv,
22592327 double**** diatomicOverlapAOs1stDeriv,
22602328 double*** tmpRotMat,
22612329 double**** tmpRotMat1stDerivs,
2262- double***** tmpDiatomicTwoElecTwoCore)const{
2263- MallocerFreer::GetInstance()->Malloc<double>(diatomicTwoElecTwoCore1stDeriv,
2330+ double***** tmpDiatomicTwoElecsTwoCores)const{
2331+ MallocerFreer::GetInstance()->Malloc<double>(diatomicTwoElecsTwoCores1stDeriv,
22642332 dxy,
22652333 dxy,
22662334 dxy,
@@ -2277,19 +2345,19 @@ void Mndo::MallocTempMatricesStaticFirstOrderFock(double****** diatomicTwoElecTw
22772345 OrbitalType_end,
22782346 OrbitalType_end,
22792347 CartesianType_end);
2280- MallocerFreer::GetInstance()->Malloc<double>(tmpDiatomicTwoElecTwoCore,
2348+ MallocerFreer::GetInstance()->Malloc<double>(tmpDiatomicTwoElecsTwoCores,
22812349 dxy,
22822350 dxy,
22832351 dxy,
22842352 dxy);
22852353 }
22862354
2287-void Mndo::FreeTempMatricesStaticFirstOrderFock(double****** diatomicTwoElecTwoCore1stDeriv,
2355+void Mndo::FreeTempMatricesStaticFirstOrderFock(double****** diatomicTwoElecsTwoCores1stDeriv,
22882356 double**** diatomicOverlapAOs1stDeriv,
22892357 double*** tmpRotMat,
22902358 double**** tmpRotMat1stDerivs,
2291- double***** tmpDiatomicTwoElecTwoCore)const{
2292- MallocerFreer::GetInstance()->Free<double>(diatomicTwoElecTwoCore1stDeriv,
2359+ double***** tmpDiatomicTwoElecsTwoCores)const{
2360+ MallocerFreer::GetInstance()->Free<double>(diatomicTwoElecsTwoCores1stDeriv,
22932361 dxy,
22942362 dxy,
22952363 dxy,
@@ -2306,7 +2374,7 @@ void Mndo::FreeTempMatricesStaticFirstOrderFock(double****** diatomicTwoElecTwoC
23062374 OrbitalType_end,
23072375 OrbitalType_end,
23082376 CartesianType_end);
2309- MallocerFreer::GetInstance()->Free<double>(tmpDiatomicTwoElecTwoCore,
2377+ MallocerFreer::GetInstance()->Free<double>(tmpDiatomicTwoElecsTwoCores,
23102378 dxy,
23112379 dxy,
23122380 dxy,
@@ -2405,7 +2473,7 @@ void Mndo::FreeTempMatricesSolveCPHF(double*** matrixCPHF,
24052473 void Mndo::CalcForceSCFElecCoreAttractionPart(double* force,
24062474 int indexAtomA,
24072475 int indexAtomB,
2408- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs) const{
2476+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
24092477 const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
24102478 int firstAOIndexA = atomA.GetFirstAOIndex();
24112479 int lastAOIndexA = atomA.GetLastAOIndex();
@@ -2418,7 +2486,7 @@ void Mndo::CalcForceSCFElecCoreAttractionPart(double* force,
24182486 indexAtomB,
24192487 mu-firstAOIndexA,
24202488 nu-firstAOIndexA,
2421- diatomicTwoElecTwoCore1stDerivs,
2489+ diatomicTwoElecsTwoCores1stDerivs,
24222490 (CartesianType)i);
24232491 }
24242492 }
@@ -2457,7 +2525,7 @@ void Mndo::CalcForceSCFOverlapAOsPart(double* force,
24572525 void Mndo::CalcForceSCFTwoElecPart(double* force,
24582526 int indexAtomA,
24592527 int indexAtomB,
2460- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs) const{
2528+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
24612529 const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
24622530 const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
24632531 int firstAOIndexA = atomA.GetFirstAOIndex();
@@ -2472,19 +2540,19 @@ void Mndo::CalcForceSCFTwoElecPart(double* force,
24722540 force[i] -= 0.5
24732541 *this->orbitalElectronPopulation[mu][nu]
24742542 *this->orbitalElectronPopulation[lambda][sigma]
2475- *diatomicTwoElecTwoCore1stDerivs[mu-firstAOIndexA]
2476- [nu-firstAOIndexA]
2477- [lambda-firstAOIndexB]
2478- [sigma-firstAOIndexB]
2479- [(CartesianType)i];
2543+ *diatomicTwoElecsTwoCores1stDerivs[mu-firstAOIndexA]
2544+ [nu-firstAOIndexA]
2545+ [lambda-firstAOIndexB]
2546+ [sigma-firstAOIndexB]
2547+ [(CartesianType)i];
24802548 force[i] += 0.25
24812549 *this->orbitalElectronPopulation[mu][lambda]
24822550 *this->orbitalElectronPopulation[nu][sigma]
2483- *diatomicTwoElecTwoCore1stDerivs[mu-firstAOIndexA]
2484- [nu-firstAOIndexA]
2485- [lambda-firstAOIndexB]
2486- [sigma-firstAOIndexB]
2487- [(CartesianType)i];
2551+ *diatomicTwoElecsTwoCores1stDerivs[mu-firstAOIndexA]
2552+ [nu-firstAOIndexA]
2553+ [lambda-firstAOIndexB]
2554+ [sigma-firstAOIndexB]
2555+ [(CartesianType)i];
24882556 }
24892557 }
24902558 }
@@ -2496,7 +2564,7 @@ void Mndo::CalcForceExcitedStaticPart(double* force,
24962564 int elecStateIndex,
24972565 int indexAtomA,
24982566 int indexAtomB,
2499- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs) const{
2567+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
25002568 const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
25012569 const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
25022570 int firstAOIndexA = atomA.GetFirstAOIndex();
@@ -2513,11 +2581,11 @@ void Mndo::CalcForceExcitedStaticPart(double* force,
25132581 -1.0*this->etaMatrixForce[elecStateIndex][mu][lambda]
25142582 *this->etaMatrixForce[elecStateIndex][nu][sigma];
25152583 force[i] += temp
2516- *diatomicTwoElecTwoCore1stDerivs[mu-firstAOIndexA]
2517- [nu-firstAOIndexA]
2518- [lambda-firstAOIndexB]
2519- [sigma-firstAOIndexB]
2520- [i];
2584+ *diatomicTwoElecsTwoCores1stDerivs[mu-firstAOIndexA]
2585+ [nu-firstAOIndexA]
2586+ [lambda-firstAOIndexB]
2587+ [sigma-firstAOIndexB]
2588+ [i];
25212589 }
25222590 }
25232591 }
@@ -2529,7 +2597,7 @@ void Mndo::CalcForceExcitedElecCoreAttractionPart(double* force,
25292597 int elecStateIndex,
25302598 int indexAtomA,
25312599 int indexAtomB,
2532- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs) const{
2600+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
25332601 const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
25342602 int firstAOIndexA = atomA.GetFirstAOIndex();
25352603 int lastAOIndexA = atomA.GetLastAOIndex();
@@ -2542,7 +2610,7 @@ void Mndo::CalcForceExcitedElecCoreAttractionPart(double* force,
25422610 indexAtomB,
25432611 mu-firstAOIndexA,
25442612 nu-firstAOIndexA,
2545- diatomicTwoElecTwoCore1stDerivs,
2613+ diatomicTwoElecsTwoCores1stDerivs,
25462614 (CartesianType)i);
25472615 }
25482616 }
@@ -2553,7 +2621,7 @@ void Mndo::CalcForceExcitedTwoElecPart(double* force,
25532621 int elecStateIndex,
25542622 int indexAtomA,
25552623 int indexAtomB,
2556- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs) const{
2624+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
25572625 const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
25582626 const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
25592627 int firstAOIndexA = atomA.GetFirstAOIndex();
@@ -2567,19 +2635,19 @@ void Mndo::CalcForceExcitedTwoElecPart(double* force,
25672635 for(int i=0; i<CartesianType_end; i++){
25682636 force[i] -= this->zMatrixForce[elecStateIndex][mu][nu]
25692637 *this->orbitalElectronPopulation[lambda][sigma]
2570- *diatomicTwoElecTwoCore1stDerivs[mu-firstAOIndexA]
2571- [nu-firstAOIndexA]
2572- [lambda-firstAOIndexB]
2573- [sigma-firstAOIndexB]
2574- [i];
2638+ *diatomicTwoElecsTwoCores1stDerivs[mu-firstAOIndexA]
2639+ [nu-firstAOIndexA]
2640+ [lambda-firstAOIndexB]
2641+ [sigma-firstAOIndexB]
2642+ [i];
25752643 force[i] += 0.50
25762644 *this->zMatrixForce[elecStateIndex][mu][lambda]
25772645 *this->orbitalElectronPopulation[nu][sigma]
2578- *diatomicTwoElecTwoCore1stDerivs[mu-firstAOIndexA]
2579- [nu-firstAOIndexA]
2580- [lambda-firstAOIndexB]
2581- [sigma-firstAOIndexB]
2582- [(CartesianType)i];
2646+ *diatomicTwoElecsTwoCores1stDerivs[mu-firstAOIndexA]
2647+ [nu-firstAOIndexA]
2648+ [lambda-firstAOIndexB]
2649+ [sigma-firstAOIndexB]
2650+ [(CartesianType)i];
25832651 }
25842652 }
25852653 }
@@ -2608,10 +2676,10 @@ void Mndo::CalcForce(const vector<int>& elecStates){
26082676 #pragma omp parallel
26092677 {
26102678 double*** diatomicOverlapAOs1stDerivs = NULL;
2611- double***** diatomicTwoElecTwoCore1stDerivs = NULL;
2612- double** tmpRotMat = NULL;
2613- double*** tmpRotMat1stDerivs = NULL;
2614- double**** tmpDiatomicTwoElecTwoCore = NULL;
2679+ double***** diatomicTwoElecsTwoCores1stDerivs = NULL;
2680+ double** tmpRotMat = NULL;
2681+ double*** tmpRotMat1stDerivs = NULL;
2682+ double**** tmpDiatomicTwoElecsTwoCores = NULL;
26152683
26162684 double** tmpDiaOverlapAOsInDiaFrame = NULL; // diatomic overlapAOs in diatomic frame
26172685 double** tmpDiaOverlapAOs1stDerivInDiaFrame = NULL; // first derivative of the diaOverlapAOs. This derivative is related to the distance between two atoms.
@@ -2622,7 +2690,7 @@ void Mndo::CalcForce(const vector<int>& elecStates){
26222690 double* tmpVectorBC = NULL; // used in dgemmm
26232691 try{
26242692 this->MallocTempMatricesCalcForce(&diatomicOverlapAOs1stDerivs,
2625- &diatomicTwoElecTwoCore1stDerivs,
2693+ &diatomicTwoElecsTwoCores1stDerivs,
26262694 &tmpDiaOverlapAOsInDiaFrame,
26272695 &tmpDiaOverlapAOs1stDerivInDiaFrame,
26282696 &tmpRotMat,
@@ -2632,7 +2700,7 @@ void Mndo::CalcForce(const vector<int>& elecStates){
26322700 &tmpRotatedDiatomicOverlapVec,
26332701 &tmpMatrixBC,
26342702 &tmpVectorBC,
2635- &tmpDiatomicTwoElecTwoCore);
2703+ &tmpDiatomicTwoElecsTwoCores);
26362704
26372705 #pragma omp for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
26382706 for(int b=0; b<this->molecule->GetNumberAtoms(); b++){
@@ -2655,11 +2723,11 @@ void Mndo::CalcForce(const vector<int>& elecStates){
26552723 atomA,
26562724 atomB);
26572725 // calc. first derivative of two elec two core interaction
2658- this->CalcDiatomicTwoElecTwoCore1stDerivatives(diatomicTwoElecTwoCore1stDerivs,
2659- tmpRotMat,
2660- tmpRotMat1stDerivs,
2661- tmpDiatomicTwoElecTwoCore,
2662- a, b);
2726+ this->CalcDiatomicTwoElecsTwoCores1stDerivatives(diatomicTwoElecsTwoCores1stDerivs,
2727+ tmpRotMat,
2728+ tmpRotMat1stDerivs,
2729+ tmpDiatomicTwoElecsTwoCores,
2730+ a, b);
26632731
26642732 // core repulsion part
26652733 double coreRepulsion[CartesianType_end] = {0.0,0.0,0.0};
@@ -2676,7 +2744,7 @@ void Mndo::CalcForce(const vector<int>& elecStates){
26762744 this->CalcForceSCFElecCoreAttractionPart(forceElecCoreAttPart,
26772745 a,
26782746 b,
2679- diatomicTwoElecTwoCore1stDerivs);
2747+ diatomicTwoElecsTwoCores1stDerivs);
26802748 // overlapAOs part (ground state)
26812749 double forceOverlapAOsPart[CartesianType_end] = {0.0,0.0,0.0};
26822750 this->CalcForceSCFOverlapAOsPart(forceOverlapAOsPart,
@@ -2688,7 +2756,7 @@ void Mndo::CalcForce(const vector<int>& elecStates){
26882756 this->CalcForceSCFTwoElecPart(forceTwoElecPart,
26892757 a,
26902758 b,
2691- diatomicTwoElecTwoCore1stDerivs);
2759+ diatomicTwoElecsTwoCores1stDerivs);
26922760 // sum up contributions from each part (ground state)
26932761 #pragma omp critical
26942762 {
@@ -2713,7 +2781,7 @@ void Mndo::CalcForce(const vector<int>& elecStates){
27132781 n,
27142782 a,
27152783 b,
2716- diatomicTwoElecTwoCore1stDerivs);
2784+ diatomicTwoElecsTwoCores1stDerivs);
27172785 // sum up contributions from static part (excited state)
27182786 #pragma omp critical
27192787 {
@@ -2731,7 +2799,7 @@ void Mndo::CalcForce(const vector<int>& elecStates){
27312799 n,
27322800 a,
27332801 b,
2734- diatomicTwoElecTwoCore1stDerivs);
2802+ diatomicTwoElecsTwoCores1stDerivs);
27352803 // overlapAOs part (excited states)
27362804 double forceExcitedOverlapAOsPart[CartesianType_end] = {0.0,0.0,0.0};
27372805 this->CalcForceExcitedOverlapAOsPart(forceExcitedOverlapAOsPart,
@@ -2745,7 +2813,7 @@ void Mndo::CalcForce(const vector<int>& elecStates){
27452813 n,
27462814 a,
27472815 b,
2748- diatomicTwoElecTwoCore1stDerivs);
2816+ diatomicTwoElecsTwoCores1stDerivs);
27492817 // sum up contributions from response part (excited state)
27502818 #pragma omp critical
27512819 {
@@ -2767,7 +2835,7 @@ void Mndo::CalcForce(const vector<int>& elecStates){
27672835 ex.Serialize(ompErrors);
27682836 }
27692837 this->FreeTempMatricesCalcForce(&diatomicOverlapAOs1stDerivs,
2770- &diatomicTwoElecTwoCore1stDerivs,
2838+ &diatomicTwoElecsTwoCores1stDerivs,
27712839 &tmpDiaOverlapAOsInDiaFrame,
27722840 &tmpDiaOverlapAOs1stDerivInDiaFrame,
27732841 &tmpRotMat,
@@ -2777,7 +2845,7 @@ void Mndo::CalcForce(const vector<int>& elecStates){
27772845 &tmpRotatedDiatomicOverlapVec,
27782846 &tmpMatrixBC,
27792847 &tmpVectorBC,
2780- &tmpDiatomicTwoElecTwoCore);
2848+ &tmpDiatomicTwoElecsTwoCores);
27812849 } // end of omp-parallelized region
27822850 // Exception throwing for omp-region
27832851 if(!ompErrors.str().empty()){
@@ -2791,7 +2859,7 @@ void Mndo::CalcForce(const vector<int>& elecStates){
27912859 }
27922860
27932861 void Mndo::MallocTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
2794- double****** diatomicTwoElecTwoCore1stDerivs,
2862+ double****** diatomicTwoElecsTwoCores1stDerivs,
27952863 double*** tmpDiaOverlapAOsInDiaFrame,
27962864 double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
27972865 double*** tmpRotMat,
@@ -2801,12 +2869,12 @@ void Mndo::MallocTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
28012869 double** tmpRotatedDiatomicOverlapVec,
28022870 double*** tmpMatrixBC,
28032871 double** tmpVectorBC,
2804- double***** tmpDiatomicTwoElecTwoCore) const{
2872+ double***** tmpDiatomicTwoElecsTwoCores) const{
28052873 MallocerFreer::GetInstance()->Malloc<double>(diatomicOverlapAOs1stDerivs,
28062874 OrbitalType_end,
28072875 OrbitalType_end,
28082876 CartesianType_end);
2809- MallocerFreer::GetInstance()->Malloc<double>(diatomicTwoElecTwoCore1stDerivs,
2877+ MallocerFreer::GetInstance()->Malloc<double>(diatomicTwoElecsTwoCores1stDerivs,
28102878 dxy,
28112879 dxy,
28122880 dxy,
@@ -2838,7 +2906,7 @@ void Mndo::MallocTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
28382906 OrbitalType_end);
28392907 MallocerFreer::GetInstance()->Malloc<double>(tmpVectorBC,
28402908 OrbitalType_end*OrbitalType_end);
2841- MallocerFreer::GetInstance()->Malloc<double>(tmpDiatomicTwoElecTwoCore,
2909+ MallocerFreer::GetInstance()->Malloc<double>(tmpDiatomicTwoElecsTwoCores,
28422910 dxy,
28432911 dxy,
28442912 dxy,
@@ -2846,7 +2914,7 @@ void Mndo::MallocTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
28462914 }
28472915
28482916 void Mndo::FreeTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
2849- double****** diatomicTwoElecTwoCore1stDerivs,
2917+ double****** diatomicTwoElecsTwoCores1stDerivs,
28502918 double*** tmpDiaOverlapAOsInDiaFrame,
28512919 double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
28522920 double*** tmpRotMat,
@@ -2856,12 +2924,12 @@ void Mndo::FreeTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
28562924 double** tmpRotatedDiatomicOverlapVec,
28572925 double*** tmpMatrixBC,
28582926 double** tmpVectorBC,
2859- double***** tmpDiatomicTwoElecTwoCore) const{
2927+ double***** tmpDiatomicTwoElecsTwoCores) const{
28602928 MallocerFreer::GetInstance()->Free<double>(diatomicOverlapAOs1stDerivs,
28612929 OrbitalType_end,
28622930 OrbitalType_end,
28632931 CartesianType_end);
2864- MallocerFreer::GetInstance()->Free<double>(diatomicTwoElecTwoCore1stDerivs,
2932+ MallocerFreer::GetInstance()->Free<double>(diatomicTwoElecsTwoCores1stDerivs,
28652933 dxy,
28662934 dxy,
28672935 dxy,
@@ -2893,7 +2961,7 @@ void Mndo::FreeTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
28932961 OrbitalType_end);
28942962 MallocerFreer::GetInstance()->Free<double>(tmpVectorBC,
28952963 OrbitalType_end*OrbitalType_end);
2896- MallocerFreer::GetInstance()->Free<double>(tmpDiatomicTwoElecTwoCore,
2964+ MallocerFreer::GetInstance()->Free<double>(tmpDiatomicTwoElecsTwoCores,
28972965 dxy,
28982966 dxy,
28992967 dxy,
@@ -2926,12 +2994,12 @@ double Mndo::GetSmallQElement(int moI,
29262994 for(int lambda=firstAOIndexB; lambda<=lastAOIndexB; lambda++){
29272995 for(int sigma=lambda; sigma<=lastAOIndexB; sigma++){
29282996 double twoElecInt = 0.0;
2929- twoElecInt = this->twoElecTwoCore[A]
2930- [B]
2931- [mu-firstAOIndexA]
2932- [nu-firstAOIndexA]
2933- [lambda-firstAOIndexB]
2934- [sigma-firstAOIndexB];
2997+ twoElecInt = this->twoElecsTwoAtomCores[A]
2998+ [B]
2999+ [mu-firstAOIndexA]
3000+ [nu-firstAOIndexA]
3001+ [lambda-firstAOIndexB]
3002+ [sigma-firstAOIndexB];
29353003 double temp = 0.0;
29363004 if(isMoPOcc){
29373005 int p = numberOcc - (moP+1);
@@ -3201,7 +3269,7 @@ double Mndo::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
32013269 }
32023270 double gamma = 0.0;
32033271 if(A!=B){
3204- gamma = this->twoElecTwoCore[A][B][muOffSet][nuOffSet][lambdaOffSet][sigmaOffSet];
3272+ gamma = this->twoElecsTwoAtomCores[A][B][muOffSet][nuOffSet][lambdaOffSet][sigmaOffSet];
32053273 }
32063274 else{
32073275 if(mu==nu && lambda==sigma){
@@ -3289,12 +3357,12 @@ double Mndo::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
32893357 for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
32903358 twoElec[A*dxy*dxy+(mu-firstAOIndexA)*dxy+(nu-firstAOIndexA)]
32913359 [B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)] =
3292- this->twoElecTwoCore[A]
3293- [B]
3294- [mu-firstAOIndexA]
3295- [nu-firstAOIndexA]
3296- [lambda-firstAOIndexB]
3297- [sigma-firstAOIndexB];
3360+ this->twoElecsTwoAtomCores[A]
3361+ [B]
3362+ [mu-firstAOIndexA]
3363+ [nu-firstAOIndexA]
3364+ [lambda-firstAOIndexB]
3365+ [sigma-firstAOIndexB];
32983366 }
32993367 }
33003368 }
@@ -3411,12 +3479,12 @@ double Mndo::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
34113479 for(int sigma=firstAOIndexB; sigma<=lastAOIndexB; sigma++){
34123480 twoElec[A*dxy*dxy+(mu-firstAOIndexA)*dxy+(nu-firstAOIndexA)]
34133481 [B*dxy*dxy+(lambda-firstAOIndexB)*dxy+(sigma-firstAOIndexB)] =
3414- this->twoElecTwoCore[A]
3415- [B]
3416- [mu-firstAOIndexA]
3417- [nu-firstAOIndexA]
3418- [lambda-firstAOIndexB]
3419- [sigma-firstAOIndexB];
3482+ this->twoElecsTwoAtomCores[A]
3483+ [B]
3484+ [mu-firstAOIndexA]
3485+ [nu-firstAOIndexA]
3486+ [lambda-firstAOIndexB]
3487+ [sigma-firstAOIndexB];
34203488 }
34213489 }
34223490 }
@@ -3486,15 +3554,15 @@ double Mndo::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) const{
34863554 return value;
34873555 }
34883556
3489-void Mndo::CalcTwoElecTwoCore(double****** twoElecTwoCore,
3557+void Mndo::CalcTwoElecsTwoCores(double****** twoElecsTwoAtomCores,
34903558 const Molecule& molecule) const{
34913559 #ifdef MOLDS_DBG
3492- if(twoElecTwoCore == NULL){
3493- throw MolDSException(this->errorMessageCalcTwoElecTwoCoreNullMatrix);
3560+ if(twoElecsTwoAtomCores == NULL){
3561+ throw MolDSException(this->errorMessageCalcTwoElecsTwoCoresNullMatrix);
34943562 }
34953563 #endif
34963564 int totalNumberAtoms = molecule.GetNumberAtoms();
3497- MallocerFreer::GetInstance()->Initialize<double>(twoElecTwoCore,
3565+ MallocerFreer::GetInstance()->Initialize<double>(twoElecsTwoAtomCores,
34983566 totalNumberAtoms,
34993567 totalNumberAtoms,
35003568 dxy, dxy, dxy, dxy);
@@ -3512,34 +3580,34 @@ void Mndo::CalcTwoElecTwoCore(double****** twoElecTwoCore,
35123580 if(mpiRank == calcRank){
35133581 #pragma omp parallel
35143582 {
3515- double**** diatomicTwoElecTwoCore = NULL;
3516- double* tmpDiatomicTwoElecTwoCore = NULL;
3517- double** tmpRotMat = NULL;
3518- double** tmpMatrixBC = NULL;
3519- double* tmpVectorBC = NULL;
3583+ double**** diatomicTwoElecsTwoCores = NULL;
3584+ double* tmpDiatomicTwoElecsTwoCores = NULL;
3585+ double** tmpRotMat = NULL;
3586+ double** tmpMatrixBC = NULL;
3587+ double* tmpVectorBC = NULL;
35203588 try{
3521- MallocerFreer::GetInstance()->Malloc<double>(&diatomicTwoElecTwoCore, dxy, dxy, dxy, dxy);
3522- MallocerFreer::GetInstance()->Malloc<double>(&tmpDiatomicTwoElecTwoCore, dxy*dxy*dxy*dxy);
3523- MallocerFreer::GetInstance()->Malloc<double>(&tmpRotMat, OrbitalType_end, OrbitalType_end);
3524- MallocerFreer::GetInstance()->Malloc<double>(&tmpMatrixBC, dxy*dxy, dxy*dxy);
3525- MallocerFreer::GetInstance()->Malloc<double>(&tmpVectorBC, dxy*dxy*dxy*dxy);
3589+ MallocerFreer::GetInstance()->Malloc<double>(&diatomicTwoElecsTwoCores, dxy, dxy, dxy, dxy);
3590+ MallocerFreer::GetInstance()->Malloc<double>(&tmpDiatomicTwoElecsTwoCores, dxy*dxy*dxy*dxy);
3591+ MallocerFreer::GetInstance()->Malloc<double>(&tmpRotMat, OrbitalType_end, OrbitalType_end);
3592+ MallocerFreer::GetInstance()->Malloc<double>(&tmpMatrixBC, dxy*dxy, dxy*dxy);
3593+ MallocerFreer::GetInstance()->Malloc<double>(&tmpVectorBC, dxy*dxy*dxy*dxy);
35263594 // note that terms with condition a==b are not needed to calculate.
35273595 #pragma omp for schedule(dynamic, MOLDS_OMP_DYNAMIC_CHUNK_SIZE)
35283596 for(int b=a+1; b<totalNumberAtoms; b++){
3529- this->CalcDiatomicTwoElecTwoCore(diatomicTwoElecTwoCore,
3530- tmpDiatomicTwoElecTwoCore,
3531- tmpRotMat,
3532- tmpMatrixBC,
3533- tmpVectorBC,
3534- a, b);
3597+ this->CalcDiatomicTwoElecsTwoCores(diatomicTwoElecsTwoCores,
3598+ tmpDiatomicTwoElecsTwoCores,
3599+ tmpRotMat,
3600+ tmpMatrixBC,
3601+ tmpVectorBC,
3602+ a, b);
35353603 int i=0;
35363604 for(int mu=0; mu<dxy; mu++){
35373605 for(int nu=mu; nu<dxy; nu++){
35383606 int j=0;
35393607 for(int lambda=0; lambda<dxy; lambda++){
35403608 for(int sigma=lambda; sigma<dxy; sigma++){
3541- this->twoElecTwoCoreMpiBuff[a][b][i][j]
3542- = diatomicTwoElecTwoCore[mu][nu][lambda][sigma];
3609+ this->twoElecsTwoAtomCoresMpiBuff[a][b][i][j]
3610+ = diatomicTwoElecsTwoCores[mu][nu][lambda][sigma];
35433611 j++;
35443612 }
35453613 }
@@ -3552,11 +3620,11 @@ void Mndo::CalcTwoElecTwoCore(double****** twoElecTwoCore,
35523620 #pragma omp critical
35533621 ex.Serialize(errorStream);
35543622 }
3555- MallocerFreer::GetInstance()->Free<double>(&diatomicTwoElecTwoCore, dxy, dxy, dxy, dxy);
3556- MallocerFreer::GetInstance()->Free<double>(&tmpDiatomicTwoElecTwoCore, dxy*dxy*dxy*dxy);
3557- MallocerFreer::GetInstance()->Free<double>(&tmpRotMat, OrbitalType_end, OrbitalType_end);
3558- MallocerFreer::GetInstance()->Free<double>(&tmpMatrixBC, dxy*dxy, dxy*dxy);
3559- MallocerFreer::GetInstance()->Free<double>(&tmpVectorBC, dxy*dxy*dxy*dxy);
3623+ MallocerFreer::GetInstance()->Free<double>(&diatomicTwoElecsTwoCores, dxy, dxy, dxy, dxy);
3624+ MallocerFreer::GetInstance()->Free<double>(&tmpDiatomicTwoElecsTwoCores, dxy*dxy*dxy*dxy);
3625+ MallocerFreer::GetInstance()->Free<double>(&tmpRotMat, OrbitalType_end, OrbitalType_end);
3626+ MallocerFreer::GetInstance()->Free<double>(&tmpMatrixBC, dxy*dxy, dxy*dxy);
3627+ MallocerFreer::GetInstance()->Free<double>(&tmpVectorBC, dxy*dxy*dxy*dxy);
35603628 }
35613629 }
35623630 if(errorStream.str().empty()){
@@ -3565,7 +3633,7 @@ void Mndo::CalcTwoElecTwoCore(double****** twoElecTwoCore,
35653633 OrbitalType twoElecLimit = dxy;
35663634 int numBuff = (twoElecLimit+1)*twoElecLimit/2;
35673635 int num = (totalNumberAtoms-b)*numBuff*numBuff;
3568- asyncCommunicator.SetBroadcastedMessage(&this->twoElecTwoCoreMpiBuff[a][b][0][0], num, calcRank);
3636+ asyncCommunicator.SetBroadcastedMessage(&this->twoElecsTwoAtomCoresMpiBuff[a][b][0][0], num, calcRank);
35693637 }
35703638 }
35713639 }
@@ -3584,15 +3652,15 @@ void Mndo::CalcTwoElecTwoCore(double****** twoElecTwoCore,
35843652 int j=0;
35853653 for(int lambda=0; lambda<dxy; lambda++){
35863654 for(int sigma=lambda; sigma<dxy; sigma++){
3587- double value = this->twoElecTwoCoreMpiBuff[a][b][i][j];
3588- twoElecTwoCore[a][b][mu][nu][lambda][sigma] = value;
3589- twoElecTwoCore[a][b][mu][nu][sigma][lambda] = value;
3590- twoElecTwoCore[a][b][nu][mu][lambda][sigma] = value;
3591- twoElecTwoCore[a][b][nu][mu][sigma][lambda] = value;
3592- twoElecTwoCore[b][a][lambda][sigma][mu][nu] = value;
3593- twoElecTwoCore[b][a][lambda][sigma][nu][mu] = value;
3594- twoElecTwoCore[b][a][sigma][lambda][mu][nu] = value;
3595- twoElecTwoCore[b][a][sigma][lambda][nu][mu] = value;
3655+ double value = this->twoElecsTwoAtomCoresMpiBuff[a][b][i][j];
3656+ twoElecsTwoAtomCores[a][b][mu][nu][lambda][sigma] = value;
3657+ twoElecsTwoAtomCores[a][b][mu][nu][sigma][lambda] = value;
3658+ twoElecsTwoAtomCores[a][b][nu][mu][lambda][sigma] = value;
3659+ twoElecsTwoAtomCores[a][b][nu][mu][sigma][lambda] = value;
3660+ twoElecsTwoAtomCores[b][a][lambda][sigma][mu][nu] = value;
3661+ twoElecsTwoAtomCores[b][a][lambda][sigma][nu][mu] = value;
3662+ twoElecsTwoAtomCores[b][a][sigma][lambda][mu][nu] = value;
3663+ twoElecsTwoAtomCores[b][a][sigma][lambda][nu][mu] = value;
35963664 j++;
35973665 }
35983666 }
@@ -3611,18 +3679,75 @@ void Mndo::CalcTwoElecTwoCore(double****** twoElecTwoCore,
36113679 // Note that atomA != atomB.
36123680 // Note taht d-orbital cannot be treated,
36133681 // that is, matrix[dxy][dxy][dxy][dxy] cannot be treatable.
3614-void Mndo::CalcDiatomicTwoElecTwoCore(double**** matrix,
3682+void Mndo::CalcDiatomicTwoElecsTwoCoresPointCharge(double**** matrix,
36153683 double* tmpVec,
36163684 double** tmpRotMat,
36173685 double** tmpMatrixBC,
36183686 double* tmpVectorBC,
3619- int indexAtomA,
3620- int indexAtomB) const{
3687+ const Atom& atomA,
3688+ const Atom& pc) const{
3689+ MallocerFreer::GetInstance()->Initialize<double>(matrix, dxy, dxy, dxy, dxy);
3690+
3691+ // calclation in diatomic frame
3692+ for(int mu=0; mu<atomA.GetValenceSize(); mu++){
3693+ for(int nu=mu; nu<atomA.GetValenceSize(); nu++){
3694+ for(int lambda=0; lambda<pc.GetValenceSize(); lambda++){
3695+ for(int sigma=lambda; sigma<pc.GetValenceSize(); sigma++){
3696+ double value = this->GetNddoRepulsionIntegralPointCharge(
3697+ atomA,
3698+ atomA.GetValence(mu),
3699+ atomA.GetValence(nu),
3700+ pc,
3701+ pc.GetValence(lambda),
3702+ pc.GetValence(sigma));
3703+ matrix[mu][nu][lambda][sigma] = value;
3704+ matrix[mu][nu][sigma][lambda] = value;
3705+ matrix[nu][mu][lambda][sigma] = value;
3706+ matrix[nu][mu][sigma][lambda] = value;
3707+ }
3708+ }
3709+ }
3710+ }
3711+ // rotate matirix into the space frame
3712+ this->CalcRotatingMatrix(tmpRotMat, atomA, pc);
3713+ this->RotateDiatomicTwoElecsTwoCoresToSpaceFrame(matrix, tmpVec, tmpRotMat, tmpMatrixBC, tmpVectorBC);
3714+
3715+ /*
3716+ this->OutputLog("(mu, nu | lambda, sigma) matrix\n");
3717+ for(int mu=0; mu<dxy; mu++){
3718+ for(int nu=0; nu<dxy; nu++){
3719+ for(int lambda=0; lambda<dxy; lambda++){
3720+ for(int sigma=0; sigma<dxy; sigma++){
3721+ this->OutputLog(boost::format("mu=%d nu=%d lambda=%d sigma=%d $e\n") % mu
3722+ % nu
3723+ % lambda
3724+ % sigma
3725+ % matrix[mu][nu][lambda][sigma]);
3726+ }
3727+ }
3728+ }
3729+ }
3730+ */
3731+}
3732+// Calculation of two electrons two cores integral (mu, nu | lambda, sigma) in space fixed frame,
3733+// taht is, Eq. (9) in ref. [DT_1977-2].
3734+// mu and nu are included in atomA's AOs.
3735+// lambda and sigma are included in atomB's AOs.
3736+// Note that atomA != atomB.
3737+// Note taht d-orbital cannot be treated,
3738+// that is, matrix[dxy][dxy][dxy][dxy] cannot be treatable.
3739+void Mndo::CalcDiatomicTwoElecsTwoCores(double**** matrix,
3740+ double* tmpVec,
3741+ double** tmpRotMat,
3742+ double** tmpMatrixBC,
3743+ double* tmpVectorBC,
3744+ int indexAtomA,
3745+ int indexAtomB) const{
36213746 const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
36223747 const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
36233748 if(indexAtomA == indexAtomB){
36243749 stringstream ss;
3625- ss << this->errorMessageCalcDiatomicTwoElecTwoCoreSameAtoms;
3750+ ss << this->errorMessageCalcDiatomicTwoElecsTwoCoresSameAtoms;
36263751 ss << this->errorMessageAtomA << indexAtomA
36273752 << AtomTypeStr(atomA.GetAtomType()) << endl;
36283753 ss << this->errorMessageAtomB << indexAtomB
@@ -3632,7 +3757,7 @@ void Mndo::CalcDiatomicTwoElecTwoCore(double**** matrix,
36323757
36333758 #ifdef MOLDS_DBG
36343759 if(matrix == NULL){
3635- throw MolDSException(this->errorMessageCalcDiatomicTwoElecTwoCoreNullMatrix);
3760+ throw MolDSException(this->errorMessageCalcDiatomicTwoElecsTwoCoresNullMatrix);
36363761 }
36373762 #endif
36383763 MallocerFreer::GetInstance()->Initialize<double>(matrix, dxy, dxy, dxy, dxy);
@@ -3659,7 +3784,7 @@ void Mndo::CalcDiatomicTwoElecTwoCore(double**** matrix,
36593784 }
36603785 // rotate matirix into the space frame
36613786 this->CalcRotatingMatrix(tmpRotMat, atomA, atomB);
3662- this->RotateDiatomicTwoElecTwoCoreToSpaceFrame(matrix, tmpVec, tmpRotMat, tmpMatrixBC, tmpVectorBC);
3787+ this->RotateDiatomicTwoElecsTwoCoresToSpaceFrame(matrix, tmpVec, tmpRotMat, tmpMatrixBC, tmpVectorBC);
36633788
36643789 /*
36653790 this->OutputLog("(mu, nu | lambda, sigma) matrix\n");
@@ -3687,17 +3812,17 @@ void Mndo::CalcDiatomicTwoElecTwoCore(double**** matrix,
36873812 // Note that atomA != atomB.
36883813 // Note taht d-orbital cannot be treated,
36893814 // that is, matrix[dxy][dxy][dxy][dxy][CartesianType_end] cannot be treatable.
3690-void Mndo::CalcDiatomicTwoElecTwoCore1stDerivatives(double***** matrix,
3691- double** tmpRotMat,
3692- double*** tmpRotMat1stDerivs,
3693- double**** tmpDiatomicTwoElecTwoCore,
3694- int indexAtomA,
3695- int indexAtomB) const{
3815+void Mndo::CalcDiatomicTwoElecsTwoCores1stDerivatives(double***** matrix,
3816+ double** tmpRotMat,
3817+ double*** tmpRotMat1stDerivs,
3818+ double**** tmpDiatomicTwoElecsTwoCores,
3819+ int indexAtomA,
3820+ int indexAtomB) const{
36963821 const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
36973822 const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
36983823 if(indexAtomA == indexAtomB){
36993824 stringstream ss;
3700- ss << this->errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesSameAtoms;
3825+ ss << this->errorMessageCalcDiatomicTwoElecsTwoCores1stDerivativesSameAtoms;
37013826 ss << this->errorMessageAtomA << indexAtomA
37023827 << AtomTypeStr(atomA.GetAtomType()) << endl;
37033828 ss << this->errorMessageAtomB << indexAtomB
@@ -3707,7 +3832,7 @@ void Mndo::CalcDiatomicTwoElecTwoCore1stDerivatives(double***** matrix,
37073832
37083833 #ifdef MOLDS_DBG
37093834 if(matrix == NULL){
3710- throw MolDSException(this->errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesNullMatrix);
3835+ throw MolDSException(this->errorMessageCalcDiatomicTwoElecsTwoCores1stDerivativesNullMatrix);
37113836 }
37123837 #endif
37133838 MallocerFreer::GetInstance()->Initialize<double>(matrix,
@@ -3736,7 +3861,7 @@ void Mndo::CalcDiatomicTwoElecTwoCore1stDerivatives(double***** matrix,
37363861 matrix[nu][mu][sigma][lambda][dimA] = matrix[mu][nu][lambda][sigma][dimA];
37373862 matrix[mu][nu][sigma][lambda][dimA] = matrix[mu][nu][lambda][sigma][dimA];
37383863 }
3739- tmpDiatomicTwoElecTwoCore[mu][nu][lambda][sigma]
3864+ tmpDiatomicTwoElecsTwoCores[mu][nu][lambda][sigma]
37403865 = this->GetNddoRepulsionIntegral(
37413866 atomA,
37423867 atomA.GetValence(mu),
@@ -3744,9 +3869,9 @@ void Mndo::CalcDiatomicTwoElecTwoCore1stDerivatives(double***** matrix,
37443869 atomB,
37453870 atomB.GetValence(lambda),
37463871 atomB.GetValence(sigma));
3747- tmpDiatomicTwoElecTwoCore[nu][mu][lambda][sigma] = tmpDiatomicTwoElecTwoCore[mu][nu][lambda][sigma];
3748- tmpDiatomicTwoElecTwoCore[nu][mu][sigma][lambda] = tmpDiatomicTwoElecTwoCore[mu][nu][lambda][sigma];
3749- tmpDiatomicTwoElecTwoCore[mu][nu][sigma][lambda] = tmpDiatomicTwoElecTwoCore[mu][nu][lambda][sigma];
3872+ tmpDiatomicTwoElecsTwoCores[nu][mu][lambda][sigma] = tmpDiatomicTwoElecsTwoCores[mu][nu][lambda][sigma];
3873+ tmpDiatomicTwoElecsTwoCores[nu][mu][sigma][lambda] = tmpDiatomicTwoElecsTwoCores[mu][nu][lambda][sigma];
3874+ tmpDiatomicTwoElecsTwoCores[mu][nu][sigma][lambda] = tmpDiatomicTwoElecsTwoCores[mu][nu][lambda][sigma];
37503875 }
37513876 }
37523877 }
@@ -3755,10 +3880,10 @@ void Mndo::CalcDiatomicTwoElecTwoCore1stDerivatives(double***** matrix,
37553880 // rotate matirix into the space frame
37563881 this->CalcRotatingMatrix(tmpRotMat, atomA, atomB);
37573882 this->CalcRotatingMatrix1stDerivatives(tmpRotMat1stDerivs, atomA, atomB);
3758- this->RotateDiatomicTwoElecTwoCore1stDerivativesToSpaceFrame(matrix,
3759- tmpDiatomicTwoElecTwoCore,
3760- tmpRotMat,
3761- tmpRotMat1stDerivs);
3883+ this->RotateDiatomicTwoElecsTwoCores1stDerivativesToSpaceFrame(matrix,
3884+ tmpDiatomicTwoElecsTwoCores,
3885+ tmpRotMat,
3886+ tmpRotMat1stDerivs);
37623887 }
37633888
37643889 // Calculation of second derivatives of the two electrons two cores integral in space fixed frame,
@@ -3769,19 +3894,19 @@ void Mndo::CalcDiatomicTwoElecTwoCore1stDerivatives(double***** matrix,
37693894 // Note that atomA != atomB.
37703895 // Note taht d-orbital cannot be treated,
37713896 // that is, matrix[dxy][dxy][dxy][dxy][CartesianType_end][CartesianType_end] cannot be treatable.
3772-void Mndo::CalcDiatomicTwoElecTwoCore2ndDerivatives(double****** matrix,
3773- double** tmpRotMat,
3774- double*** tmpRotMat1stDerivs,
3775- double**** tmpRotMat2ndDerivs,
3776- double**** tmpDiatomicTwoElecTwoCore,
3777- double***** tmpDiatomicTwoElecTwoCore1stDerivs,
3778- int indexAtomA,
3779- int indexAtomB) const{
3897+void Mndo::CalcDiatomicTwoElecsTwoCores2ndDerivatives(double****** matrix,
3898+ double** tmpRotMat,
3899+ double*** tmpRotMat1stDerivs,
3900+ double**** tmpRotMat2ndDerivs,
3901+ double**** tmpDiatomicTwoElecsTwoCores,
3902+ double***** tmpDiatomicTwoElecsTwoCores1stDerivs,
3903+ int indexAtomA,
3904+ int indexAtomB) const{
37803905 const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
37813906 const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
37823907 if(indexAtomA == indexAtomB){
37833908 stringstream ss;
3784- ss << this->errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesSameAtoms;
3909+ ss << this->errorMessageCalcDiatomicTwoElecsTwoCores2ndDerivativesSameAtoms;
37853910 ss << this->errorMessageAtomA << indexAtomA
37863911 << AtomTypeStr(atomA.GetAtomType()) << endl;
37873912 ss << this->errorMessageAtomB << indexAtomB
@@ -3791,7 +3916,7 @@ void Mndo::CalcDiatomicTwoElecTwoCore2ndDerivatives(double****** matrix,
37913916
37923917 #ifdef MOLDS_DBG
37933918 if(matrix == NULL){
3794- throw MolDSException(this->errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesNullMatrix);
3919+ throw MolDSException(this->errorMessageCalcDiatomicTwoElecsTwoCores2ndDerivativesNullMatrix);
37953920 }
37963921 #endif
37973922 MallocerFreer::GetInstance()->Initialize<double>(matrix,
@@ -3820,7 +3945,7 @@ void Mndo::CalcDiatomicTwoElecTwoCore2ndDerivatives(double****** matrix,
38203945 static_cast<CartesianType>(dimA1),
38213946 static_cast<CartesianType>(dimA2));
38223947 }
3823- tmpDiatomicTwoElecTwoCore1stDerivs[mu][nu][lambda][sigma][dimA1]
3948+ tmpDiatomicTwoElecsTwoCores1stDerivs[mu][nu][lambda][sigma][dimA1]
38243949 = this->GetNddoRepulsionIntegral1stDerivative(
38253950 atomA,
38263951 atomA.GetValence(mu),
@@ -3830,7 +3955,7 @@ void Mndo::CalcDiatomicTwoElecTwoCore2ndDerivatives(double****** matrix,
38303955 atomB.GetValence(sigma),
38313956 static_cast<CartesianType>(dimA1));
38323957 }
3833- tmpDiatomicTwoElecTwoCore[mu][nu][lambda][sigma]
3958+ tmpDiatomicTwoElecsTwoCores[mu][nu][lambda][sigma]
38343959 = this->GetNddoRepulsionIntegral(
38353960 atomA,
38363961 atomA.GetValence(mu),
@@ -3847,21 +3972,21 @@ void Mndo::CalcDiatomicTwoElecTwoCore2ndDerivatives(double****** matrix,
38473972 this->CalcRotatingMatrix(tmpRotMat, atomA, atomB);
38483973 this->CalcRotatingMatrix1stDerivatives(tmpRotMat1stDerivs, atomA, atomB);
38493974 this->CalcRotatingMatrix2ndDerivatives(tmpRotMat2ndDerivs, atomA, atomB);
3850- this->RotateDiatomicTwoElecTwoCore2ndDerivativesToSpaceFrame(matrix,
3851- tmpDiatomicTwoElecTwoCore,
3852- tmpDiatomicTwoElecTwoCore1stDerivs,
3853- tmpRotMat,
3854- tmpRotMat1stDerivs,
3855- tmpRotMat2ndDerivs);
3975+ this->RotateDiatomicTwoElecsTwoCores2ndDerivativesToSpaceFrame(matrix,
3976+ tmpDiatomicTwoElecsTwoCores,
3977+ tmpDiatomicTwoElecsTwoCores1stDerivs,
3978+ tmpRotMat,
3979+ tmpRotMat1stDerivs,
3980+ tmpRotMat2ndDerivs);
38563981 }
38573982
38583983 // Rotate 4-dimensional matrix from diatomic frame to space frame
38593984 // Note tha in this method d-orbitals can not be treatable.
3860-void Mndo::RotateDiatomicTwoElecTwoCoreToSpaceFrame(double**** matrix,
3861- double* tmpVec,
3862- double const* const* rotatingMatrix,
3863- double** tmpMatrixBC,
3864- double* tmpVectorBC) const{
3985+void Mndo::RotateDiatomicTwoElecsTwoCoresToSpaceFrame(double**** matrix,
3986+ double* tmpVec,
3987+ double const* const* rotatingMatrix,
3988+ double** tmpMatrixBC,
3989+ double* tmpVectorBC) const{
38653990 double oldMatrix[dxy][dxy][dxy][dxy];
38663991 MolDS_wrappers::Blas::GetInstance()->Dcopy(dxy*dxy*dxy*dxy, &matrix[0][0][0][0], &oldMatrix[0][0][0][0]);
38673992
@@ -3931,9 +4056,9 @@ void Mndo::RotateDiatomicTwoElecTwoCoreToSpaceFrame(double**** matrix,
39314056
39324057 // Rotate 5-dimensional matrix from diatomic frame to space frame
39334058 // Note tha in this method d-orbitals can not be treatable.
3934-void Mndo::RotateDiatomicTwoElecTwoCore1stDerivativesToSpaceFrame(
4059+void Mndo::RotateDiatomicTwoElecsTwoCores1stDerivativesToSpaceFrame(
39354060 double***** matrix,
3936- double const* const*const* const* diatomicTwoElecTwoCore,
4061+ double const* const*const* const* diatomicTwoElecsTwoCores,
39374062 double const* const* rotatingMatrix,
39384063 double const* const* const* rotMat1stDerivatives) const{
39394064
@@ -3953,15 +4078,15 @@ void Mndo::RotateDiatomicTwoElecTwoCore1stDerivativesToSpaceFrame(
39534078 double* tmpVector = NULL;
39544079 double** ptrDiatomic = NULL;
39554080 try{
3956- this->MallocTempMatricesRotateDiatomicTwoElecTwoCore1stDerivs(&twiceRotatingMatrix,
3957- &twiceRotatingMatrixDerivA,
3958- &twiceRotatingMatrixDerivB,
3959- &oldMatrix,
3960- &rotatedMatrix,
3961- &tmpRotatedVec,
3962- &tmpMatrix,
3963- &tmpVector,
3964- &ptrDiatomic);
4081+ this->MallocTempMatricesRotateDiatomicTwoElecsTwoCores1stDerivs(&twiceRotatingMatrix,
4082+ &twiceRotatingMatrixDerivA,
4083+ &twiceRotatingMatrixDerivB,
4084+ &oldMatrix,
4085+ &rotatedMatrix,
4086+ &tmpRotatedVec,
4087+ &tmpMatrix,
4088+ &tmpVector,
4089+ &ptrDiatomic);
39654090 for(int mu=0; mu<dxy; mu++){
39664091 for(int nu=0; nu<dxy; nu++){
39674092 int i=mu*dxy+nu;
@@ -3972,7 +4097,7 @@ void Mndo::RotateDiatomicTwoElecTwoCore1stDerivativesToSpaceFrame(
39724097 *rotatingMatrix[nu][sigma ];
39734098 }
39744099 }
3975- ptrDiatomic[i] = const_cast<double*>(&diatomicTwoElecTwoCore[mu][nu][0][0]);
4100+ ptrDiatomic[i] = const_cast<double*>(&diatomicTwoElecsTwoCores[mu][nu][0][0]);
39764101 }
39774102 }
39784103 for(int axis=0; axis<CartesianType_end; axis++){
@@ -4067,26 +4192,26 @@ void Mndo::RotateDiatomicTwoElecTwoCore1stDerivativesToSpaceFrame(
40674192 }
40684193 }
40694194 catch(MolDSException ex){
4070- this->FreeTempMatricesRotateDiatomicTwoElecTwoCore1stDerivs(&twiceRotatingMatrix,
4071- &twiceRotatingMatrixDerivA,
4072- &twiceRotatingMatrixDerivB,
4073- &oldMatrix,
4074- &rotatedMatrix,
4075- &tmpRotatedVec,
4076- &tmpMatrix,
4077- &tmpVector,
4078- &ptrDiatomic);
4195+ this->FreeTempMatricesRotateDiatomicTwoElecsTwoCores1stDerivs(&twiceRotatingMatrix,
4196+ &twiceRotatingMatrixDerivA,
4197+ &twiceRotatingMatrixDerivB,
4198+ &oldMatrix,
4199+ &rotatedMatrix,
4200+ &tmpRotatedVec,
4201+ &tmpMatrix,
4202+ &tmpVector,
4203+ &ptrDiatomic);
40794204 throw ex;
40804205 }
4081- this->FreeTempMatricesRotateDiatomicTwoElecTwoCore1stDerivs(&twiceRotatingMatrix,
4082- &twiceRotatingMatrixDerivA,
4083- &twiceRotatingMatrixDerivB,
4084- &oldMatrix,
4085- &rotatedMatrix,
4086- &tmpRotatedVec,
4087- &tmpMatrix,
4088- &tmpVector,
4089- &ptrDiatomic);
4206+ this->FreeTempMatricesRotateDiatomicTwoElecsTwoCores1stDerivs(&twiceRotatingMatrix,
4207+ &twiceRotatingMatrixDerivA,
4208+ &twiceRotatingMatrixDerivB,
4209+ &oldMatrix,
4210+ &rotatedMatrix,
4211+ &tmpRotatedVec,
4212+ &tmpMatrix,
4213+ &tmpVector,
4214+ &ptrDiatomic);
40904215
40914216 /*
40924217 // rotate (slow algorithm)
@@ -4107,25 +4232,25 @@ void Mndo::RotateDiatomicTwoElecTwoCore1stDerivativesToSpaceFrame(
41074232 *rotatingMatrix[lambda][k]
41084233 *rotatingMatrix[sigma][l];
41094234 matrix[mu][nu][lambda][sigma][c]
4110- += diatomicTwoElecTwoCore[i][j][k][l]
4235+ += diatomicTwoElecsTwoCores[i][j][k][l]
41114236 *rotMat1stDerivatives[mu][i][c]
41124237 *rotatingMatrix[nu][j]
41134238 *rotatingMatrix[lambda][k]
41144239 *rotatingMatrix[sigma][l];
41154240 matrix[mu][nu][lambda][sigma][c]
4116- += diatomicTwoElecTwoCore[i][j][k][l]
4241+ += diatomicTwoElecsTwoCores[i][j][k][l]
41174242 *rotatingMatrix[mu][i]
41184243 *rotMat1stDerivatives[nu][j][c]
41194244 *rotatingMatrix[lambda][k]
41204245 *rotatingMatrix[sigma][l];
41214246 matrix[mu][nu][lambda][sigma][c]
4122- += diatomicTwoElecTwoCore[i][j][k][l]
4247+ += diatomicTwoElecsTwoCores[i][j][k][l]
41234248 *rotatingMatrix[mu][i]
41244249 *rotatingMatrix[nu][j]
41254250 *rotMat1stDerivatives[lambda][k][c]
41264251 *rotatingMatrix[sigma][l];
41274252 matrix[mu][nu][lambda][sigma][c]
4128- += diatomicTwoElecTwoCore[i][j][k][l]
4253+ += diatomicTwoElecsTwoCores[i][j][k][l]
41294254 *rotatingMatrix[mu][i]
41304255 *rotatingMatrix[nu][j]
41314256 *rotatingMatrix[lambda][k]
@@ -4142,15 +4267,15 @@ void Mndo::RotateDiatomicTwoElecTwoCore1stDerivativesToSpaceFrame(
41424267 */
41434268 }
41444269
4145-void Mndo::MallocTempMatricesRotateDiatomicTwoElecTwoCore1stDerivs(double*** twiceRotatingMatrix,
4146- double*** twiceRotatingMatrixDerivA,
4147- double*** twiceRotatingMatrixDerivB,
4148- double*** oldMatrix,
4149- double*** rotatedMatrix,
4150- double** tmpRotatedVec,
4151- double*** tmpMatrix,
4152- double** tmpVector,
4153- double*** ptrDiatomic) const{
4270+void Mndo::MallocTempMatricesRotateDiatomicTwoElecsTwoCores1stDerivs(double*** twiceRotatingMatrix,
4271+ double*** twiceRotatingMatrixDerivA,
4272+ double*** twiceRotatingMatrixDerivB,
4273+ double*** oldMatrix,
4274+ double*** rotatedMatrix,
4275+ double** tmpRotatedVec,
4276+ double*** tmpMatrix,
4277+ double** tmpVector,
4278+ double*** ptrDiatomic) const{
41544279 MallocerFreer::GetInstance()->Malloc<double>(twiceRotatingMatrix, dxy*dxy, dxy*dxy);
41554280 MallocerFreer::GetInstance()->Malloc<double>(twiceRotatingMatrixDerivA, dxy*dxy, dxy*dxy);
41564281 MallocerFreer::GetInstance()->Malloc<double>(twiceRotatingMatrixDerivB, dxy*dxy, dxy*dxy);
@@ -4162,15 +4287,15 @@ void Mndo::MallocTempMatricesRotateDiatomicTwoElecTwoCore1stDerivs(double*** twi
41624287 MallocerFreer::GetInstance()->Malloc<double*>(ptrDiatomic, dxy*dxy);
41634288 }
41644289
4165-void Mndo::FreeTempMatricesRotateDiatomicTwoElecTwoCore1stDerivs(double*** twiceRotatingMatrix,
4166- double*** twiceRotatingMatrixDerivA,
4167- double*** twiceRotatingMatrixDerivB,
4168- double*** oldMatrix,
4169- double*** rotatedMatrix,
4170- double** tmpRotatedVec,
4171- double*** tmpMatrix,
4172- double** tmpVector,
4173- double*** ptrDiatomic) const{
4290+void Mndo::FreeTempMatricesRotateDiatomicTwoElecsTwoCores1stDerivs(double*** twiceRotatingMatrix,
4291+ double*** twiceRotatingMatrixDerivA,
4292+ double*** twiceRotatingMatrixDerivB,
4293+ double*** oldMatrix,
4294+ double*** rotatedMatrix,
4295+ double** tmpRotatedVec,
4296+ double*** tmpMatrix,
4297+ double** tmpVector,
4298+ double*** ptrDiatomic) const{
41744299 MallocerFreer::GetInstance()->Free<double>(twiceRotatingMatrix, dxy*dxy, dxy*dxy);
41754300 MallocerFreer::GetInstance()->Free<double>(twiceRotatingMatrixDerivA, dxy*dxy, dxy*dxy);
41764301 MallocerFreer::GetInstance()->Free<double>(twiceRotatingMatrixDerivB, dxy*dxy, dxy*dxy);
@@ -4184,10 +4309,10 @@ void Mndo::FreeTempMatricesRotateDiatomicTwoElecTwoCore1stDerivs(double*** twice
41844309
41854310 // Rotate 6-dimensional matrix from diatomic frame to space frame
41864311 // Note tha in this method d-orbitals can not be treatable.
4187-void Mndo::RotateDiatomicTwoElecTwoCore2ndDerivativesToSpaceFrame(
4312+void Mndo::RotateDiatomicTwoElecsTwoCores2ndDerivativesToSpaceFrame(
41884313 double****** matrix,
4189- double const* const* const* const* diatomicTwoElecTwoCore,
4190- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivatives,
4314+ double const* const* const* const* diatomicTwoElecsTwoCores,
4315+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivatives,
41914316 double const* const* rotatingMatrix,
41924317 double const* const* const* rotMat1stDerivatives,
41934318 double const* const* const* const* rotMat2ndDerivatives) const{
@@ -4231,18 +4356,18 @@ void Mndo::RotateDiatomicTwoElecTwoCore2ndDerivativesToSpaceFrame(
42314356 MallocerFreer::GetInstance()->Initialize<double>(tempIJK, numberTerms);
42324357 for(int l=s; l<dxy; l++){
42334358
4234- tempIJK[0] += oldMatrix [i][j][k][l][dimA1][dimA2]*rotatingMatrix [sigma][l];
4235- tempIJK[1] += diatomicTwoElecTwoCore[i][j][k][l] *rotatingMatrix [sigma][l];
4236- tempIJK[4] += diatomicTwoElecTwoCore[i][j][k][l] *rotMat2ndDerivatives[sigma][l][dimA1][dimA2];
4359+ tempIJK[0] += oldMatrix [i][j][k][l][dimA1][dimA2]*rotatingMatrix [sigma][l];
4360+ tempIJK[1] += diatomicTwoElecsTwoCores[i][j][k][l] *rotatingMatrix [sigma][l];
4361+ tempIJK[4] += diatomicTwoElecsTwoCores[i][j][k][l] *rotMat2ndDerivatives[sigma][l][dimA1][dimA2];
42374362
4238- tempIJK[5] += diatomicTwoElecTwoCore1stDerivatives[i][j][k][l][dimA1]*rotatingMatrix [sigma][l];
4239- tempIJK[8] += diatomicTwoElecTwoCore1stDerivatives[i][j][k][l][dimA1]*rotMat1stDerivatives[sigma][l][dimA2];
4363+ tempIJK[5] += diatomicTwoElecsTwoCores1stDerivatives[i][j][k][l][dimA1]*rotatingMatrix [sigma][l];
4364+ tempIJK[8] += diatomicTwoElecsTwoCores1stDerivatives[i][j][k][l][dimA1]*rotMat1stDerivatives[sigma][l][dimA2];
42404365
4241- tempIJK[9] += diatomicTwoElecTwoCore1stDerivatives[i][j][k][l][dimA2]*rotatingMatrix [sigma][l];
4242- tempIJK[12] += diatomicTwoElecTwoCore [i][j][k][l] *rotMat1stDerivatives[sigma][l][dimA2];
4366+ tempIJK[9] += diatomicTwoElecsTwoCores1stDerivatives[i][j][k][l][dimA2]*rotatingMatrix [sigma][l];
4367+ tempIJK[12] += diatomicTwoElecsTwoCores [i][j][k][l] *rotMat1stDerivatives[sigma][l][dimA2];
42434368
4244- tempIJK[21] += diatomicTwoElecTwoCore1stDerivatives[i][j][k][l][dimA2]*rotMat1stDerivatives[sigma][l][dimA1];
4245- tempIJK[22] += diatomicTwoElecTwoCore [i][j][k][l] *rotMat1stDerivatives[sigma][l][dimA1];
4369+ tempIJK[21] += diatomicTwoElecsTwoCores1stDerivatives[i][j][k][l][dimA2]*rotMat1stDerivatives[sigma][l][dimA1];
4370+ tempIJK[22] += diatomicTwoElecsTwoCores [i][j][k][l] *rotMat1stDerivatives[sigma][l][dimA1];
42464371 }
42474372 tempIJ[0] += tempIJK[0] *rotatingMatrix [lambda][k];
42484373 tempIJ[1] += tempIJK[1] *rotatingMatrix [lambda][k];
@@ -4353,154 +4478,154 @@ void Mndo::RotateDiatomicTwoElecTwoCore2ndDerivativesToSpaceFrame(
43534478 *rotatingMatrix [lambda][k]
43544479 *rotatingMatrix [sigma ][l];
43554480 // term1
4356- value += diatomicTwoElecTwoCore[i][j][k][l]
4357- *rotMat2ndDerivatives[mu ][i][dimA1][dimA2]
4481+ value += diatomicTwoElecsTwoCores[i][j][k][l]
4482+ *rotMat2ndDerivatives [mu ][i][dimA1][dimA2]
43584483 *rotatingMatrix [nu ][j]
43594484 *rotatingMatrix [lambda][k]
43604485 *rotatingMatrix [sigma ][l];
43614486 // term2
4362- value += diatomicTwoElecTwoCore[i][j][k][l]
4487+ value += diatomicTwoElecsTwoCores[i][j][k][l]
43634488 *rotatingMatrix [mu ][i]
4364- *rotMat2ndDerivatives[nu ][j][dimA1][dimA2]
4489+ *rotMat2ndDerivatives [nu ][j][dimA1][dimA2]
43654490 *rotatingMatrix [lambda][k]
43664491 *rotatingMatrix [sigma ][l];
43674492 // term3
4368- value += diatomicTwoElecTwoCore[i][j][k][l]
4493+ value += diatomicTwoElecsTwoCores[i][j][k][l]
43694494 *rotatingMatrix [mu ][i]
43704495 *rotatingMatrix [nu ][j]
4371- *rotMat2ndDerivatives[lambda][k][dimA1][dimA2]
4496+ *rotMat2ndDerivatives [lambda][k][dimA1][dimA2]
43724497 *rotatingMatrix [sigma ][l];
43734498 // term4
4374- value += diatomicTwoElecTwoCore[i][j][k][l]
4499+ value += diatomicTwoElecsTwoCores[i][j][k][l]
43754500 *rotatingMatrix [mu ][i]
43764501 *rotatingMatrix [nu ][j]
43774502 *rotatingMatrix [lambda][k]
4378- *rotMat2ndDerivatives[sigma ][l][dimA1][dimA2];
4503+ *rotMat2ndDerivatives [sigma ][l][dimA1][dimA2];
43794504
43804505 // term5
4381- value += diatomicTwoElecTwoCore1stDerivatives[i][j][k][l][dimA1]
4382- *rotMat1stDerivatives[mu ][i][dimA2]
4383- *rotatingMatrix [nu ][j]
4384- *rotatingMatrix [lambda][k]
4385- *rotatingMatrix [sigma ][l];
4506+ value += diatomicTwoElecsTwoCores1stDerivatives[i][j][k][l][dimA1]
4507+ *rotMat1stDerivatives [mu ][i][dimA2]
4508+ *rotatingMatrix [nu ][j]
4509+ *rotatingMatrix [lambda][k]
4510+ *rotatingMatrix [sigma ][l];
43864511 // term6
4387- value += diatomicTwoElecTwoCore1stDerivatives[i][j][k][l][dimA1]
4388- *rotatingMatrix [mu ][i]
4389- *rotMat1stDerivatives[nu ][j][dimA2]
4390- *rotatingMatrix [lambda][k]
4391- *rotatingMatrix [sigma ][l];
4512+ value += diatomicTwoElecsTwoCores1stDerivatives[i][j][k][l][dimA1]
4513+ *rotatingMatrix [mu ][i]
4514+ *rotMat1stDerivatives [nu ][j][dimA2]
4515+ *rotatingMatrix [lambda][k]
4516+ *rotatingMatrix [sigma ][l];
43924517 // term7
4393- value += diatomicTwoElecTwoCore1stDerivatives[i][j][k][l][dimA1]
4394- *rotatingMatrix [mu ][i]
4395- *rotatingMatrix [nu ][j]
4396- *rotMat1stDerivatives[lambda][k][dimA2]
4397- *rotatingMatrix [sigma ][l];
4518+ value += diatomicTwoElecsTwoCores1stDerivatives[i][j][k][l][dimA1]
4519+ *rotatingMatrix [mu ][i]
4520+ *rotatingMatrix [nu ][j]
4521+ *rotMat1stDerivatives [lambda][k][dimA2]
4522+ *rotatingMatrix [sigma ][l];
43984523 // term8
4399- value += diatomicTwoElecTwoCore1stDerivatives[i][j][k][l][dimA1]
4400- *rotatingMatrix [mu ][i]
4401- *rotatingMatrix [nu ][j]
4402- *rotatingMatrix [lambda][k]
4403- *rotMat1stDerivatives[sigma ][l][dimA2];
4524+ value += diatomicTwoElecsTwoCores1stDerivatives[i][j][k][l][dimA1]
4525+ *rotatingMatrix [mu ][i]
4526+ *rotatingMatrix [nu ][j]
4527+ *rotatingMatrix [lambda][k]
4528+ *rotMat1stDerivatives [sigma ][l][dimA2];
44044529
44054530 // term9
4406- value += diatomicTwoElecTwoCore1stDerivatives[i][j][k][l][dimA2]
4407- *rotMat1stDerivatives[mu ][i][dimA1]
4408- *rotatingMatrix [nu ][j]
4409- *rotatingMatrix [lambda][k]
4410- *rotatingMatrix [sigma ][l];
4531+ value += diatomicTwoElecsTwoCores1stDerivatives[i][j][k][l][dimA2]
4532+ *rotMat1stDerivatives [mu ][i][dimA1]
4533+ *rotatingMatrix [nu ][j]
4534+ *rotatingMatrix [lambda][k]
4535+ *rotatingMatrix [sigma ][l];
44114536 // term10
4412- value += diatomicTwoElecTwoCore[i][j][k][l]
4413- *rotMat1stDerivatives[mu ][i][dimA1]
4414- *rotMat1stDerivatives[nu ][j][dimA2]
4415- *rotatingMatrix [lambda][k]
4416- *rotatingMatrix [sigma ][l];
4537+ value += diatomicTwoElecsTwoCores[i][j][k][l]
4538+ *rotMat1stDerivatives [mu ][i][dimA1]
4539+ *rotMat1stDerivatives [nu ][j][dimA2]
4540+ *rotatingMatrix [lambda][k]
4541+ *rotatingMatrix [sigma ][l];
44174542 // term11
4418- value += diatomicTwoElecTwoCore[i][j][k][l]
4419- *rotMat1stDerivatives[mu ][i][dimA1]
4420- *rotatingMatrix [nu ][j]
4421- *rotMat1stDerivatives[lambda][k][dimA2]
4422- *rotatingMatrix [sigma ][l];
4543+ value += diatomicTwoElecsTwoCores[i][j][k][l]
4544+ *rotMat1stDerivatives [mu ][i][dimA1]
4545+ *rotatingMatrix [nu ][j]
4546+ *rotMat1stDerivatives [lambda][k][dimA2]
4547+ *rotatingMatrix [sigma ][l];
44234548 // term12
4424- value += diatomicTwoElecTwoCore[i][j][k][l]
4425- *rotMat1stDerivatives[mu ][i][dimA1]
4426- *rotatingMatrix [nu ][j]
4427- *rotatingMatrix [lambda][k]
4428- *rotMat1stDerivatives[sigma ][l][dimA2];
4549+ value += diatomicTwoElecsTwoCores[i][j][k][l]
4550+ *rotMat1stDerivatives [mu ][i][dimA1]
4551+ *rotatingMatrix [nu ][j]
4552+ *rotatingMatrix [lambda][k]
4553+ *rotMat1stDerivatives [sigma ][l][dimA2];
44294554
44304555 // term13
4431- value += diatomicTwoElecTwoCore1stDerivatives[i][j][k][l][dimA2]
4432- *rotatingMatrix [mu ][i]
4433- *rotMat1stDerivatives[nu ][j][dimA1]
4434- *rotatingMatrix [lambda][k]
4435- *rotatingMatrix [sigma ][l];
4556+ value += diatomicTwoElecsTwoCores1stDerivatives[i][j][k][l][dimA2]
4557+ *rotatingMatrix [mu ][i]
4558+ *rotMat1stDerivatives [nu ][j][dimA1]
4559+ *rotatingMatrix [lambda][k]
4560+ *rotatingMatrix [sigma ][l];
44364561 // term14
4437- value += diatomicTwoElecTwoCore[i][j][k][l]
4438- *rotMat1stDerivatives[mu ][i][dimA2]
4439- *rotMat1stDerivatives[nu ][j][dimA1]
4440- *rotatingMatrix [lambda][k]
4441- *rotatingMatrix [sigma ][l];
4562+ value += diatomicTwoElecsTwoCores[i][j][k][l]
4563+ *rotMat1stDerivatives [mu ][i][dimA2]
4564+ *rotMat1stDerivatives [nu ][j][dimA1]
4565+ *rotatingMatrix [lambda][k]
4566+ *rotatingMatrix [sigma ][l];
44424567 // term15
4443- value += diatomicTwoElecTwoCore[i][j][k][l]
4444- *rotatingMatrix [mu ][i]
4445- *rotMat1stDerivatives[nu ][j][dimA1]
4446- *rotMat1stDerivatives[lambda][k][dimA2]
4447- *rotatingMatrix [sigma ][l];
4568+ value += diatomicTwoElecsTwoCores[i][j][k][l]
4569+ *rotatingMatrix [mu ][i]
4570+ *rotMat1stDerivatives [nu ][j][dimA1]
4571+ *rotMat1stDerivatives [lambda][k][dimA2]
4572+ *rotatingMatrix [sigma ][l];
44484573 // term16
4449- value += diatomicTwoElecTwoCore[i][j][k][l]
4450- *rotatingMatrix [mu ][i]
4451- *rotMat1stDerivatives[nu ][j][dimA1]
4452- *rotatingMatrix [lambda][k]
4453- *rotMat1stDerivatives[sigma ][l][dimA2];
4574+ value += diatomicTwoElecsTwoCores[i][j][k][l]
4575+ *rotatingMatrix [mu ][i]
4576+ *rotMat1stDerivatives [nu ][j][dimA1]
4577+ *rotatingMatrix [lambda][k]
4578+ *rotMat1stDerivatives [sigma ][l][dimA2];
44544579
44554580 // term17
4456- value += diatomicTwoElecTwoCore1stDerivatives[i][j][k][l][dimA2]
4457- *rotatingMatrix [mu ][i]
4458- *rotatingMatrix [nu ][j]
4459- *rotMat1stDerivatives[lambda][k][dimA1]
4460- *rotatingMatrix [sigma ][l];
4581+ value += diatomicTwoElecsTwoCores1stDerivatives[i][j][k][l][dimA2]
4582+ *rotatingMatrix [mu ][i]
4583+ *rotatingMatrix [nu ][j]
4584+ *rotMat1stDerivatives [lambda][k][dimA1]
4585+ *rotatingMatrix [sigma ][l];
44614586 // term18
4462- value += diatomicTwoElecTwoCore[i][j][k][l]
4463- *rotMat1stDerivatives[mu ][i][dimA2]
4464- *rotatingMatrix [nu ][j]
4465- *rotMat1stDerivatives[lambda][k][dimA1]
4466- *rotatingMatrix [sigma ][l];
4587+ value += diatomicTwoElecsTwoCores[i][j][k][l]
4588+ *rotMat1stDerivatives [mu ][i][dimA2]
4589+ *rotatingMatrix [nu ][j]
4590+ *rotMat1stDerivatives [lambda][k][dimA1]
4591+ *rotatingMatrix [sigma ][l];
44674592 // term19
4468- value += diatomicTwoElecTwoCore[i][j][k][l]
4469- *rotatingMatrix [mu ][i]
4470- *rotMat1stDerivatives[nu ][j][dimA2]
4471- *rotMat1stDerivatives[lambda][k][dimA1]
4472- *rotatingMatrix [sigma ][l];
4593+ value += diatomicTwoElecsTwoCores[i][j][k][l]
4594+ *rotatingMatrix [mu ][i]
4595+ *rotMat1stDerivatives [nu ][j][dimA2]
4596+ *rotMat1stDerivatives [lambda][k][dimA1]
4597+ *rotatingMatrix [sigma ][l];
44734598 // term20
4474- value += diatomicTwoElecTwoCore[i][j][k][l]
4475- *rotatingMatrix [mu ][i]
4476- *rotatingMatrix [nu ][j]
4477- *rotMat1stDerivatives[lambda][k][dimA1]
4478- *rotMat1stDerivatives[sigma ][l][dimA2];
4599+ value += diatomicTwoElecsTwoCores[i][j][k][l]
4600+ *rotatingMatrix [mu ][i]
4601+ *rotatingMatrix [nu ][j]
4602+ *rotMat1stDerivatives [lambda][k][dimA1]
4603+ *rotMat1stDerivatives [sigma ][l][dimA2];
44794604
44804605 // term21
4481- value += diatomicTwoElecTwoCore1stDerivatives[i][j][k][l][dimA2]
4482- *rotatingMatrix [mu ][i]
4483- *rotatingMatrix [nu ][j]
4484- *rotatingMatrix [lambda][k]
4485- *rotMat1stDerivatives[sigma ][l][dimA1];
4606+ value += diatomicTwoElecsTwoCores1stDerivatives[i][j][k][l][dimA2]
4607+ *rotatingMatrix [mu ][i]
4608+ *rotatingMatrix [nu ][j]
4609+ *rotatingMatrix [lambda][k]
4610+ *rotMat1stDerivatives [sigma ][l][dimA1];
44864611 // term22
4487- value += diatomicTwoElecTwoCore[i][j][k][l]
4488- *rotMat1stDerivatives[mu ][i][dimA2]
4489- *rotatingMatrix [nu ][j]
4490- *rotatingMatrix [lambda][k]
4491- *rotMat1stDerivatives[sigma ][l][dimA1];
4612+ value += diatomicTwoElecsTwoCores[i][j][k][l]
4613+ *rotMat1stDerivatives [mu ][i][dimA2]
4614+ *rotatingMatrix [nu ][j]
4615+ *rotatingMatrix [lambda][k]
4616+ *rotMat1stDerivatives [sigma ][l][dimA1];
44924617 // term23
4493- value += diatomicTwoElecTwoCore[i][j][k][l]
4494- *rotatingMatrix [mu ][i]
4495- *rotMat1stDerivatives[nu ][j][dimA2]
4496- *rotatingMatrix [lambda][k]
4497- *rotMat1stDerivatives[sigma ][l][dimA1];
4618+ value += diatomicTwoElecsTwoCores[i][j][k][l]
4619+ *rotatingMatrix [mu ][i]
4620+ *rotMat1stDerivatives [nu ][j][dimA2]
4621+ *rotatingMatrix [lambda][k]
4622+ *rotMat1stDerivatives [sigma ][l][dimA1];
44984623 // term24
4499- value += diatomicTwoElecTwoCore[i][j][k][l]
4500- *rotatingMatrix [mu ][i]
4501- *rotatingMatrix [nu ][j]
4502- *rotMat1stDerivatives[lambda][k][dimA2]
4503- *rotMat1stDerivatives[sigma ][l][dimA1];
4624+ value += diatomicTwoElecsTwoCores[i][j][k][l]
4625+ *rotatingMatrix [mu ][i]
4626+ *rotatingMatrix [nu ][j]
4627+ *rotMat1stDerivatives [lambda][k][dimA2]
4628+ *rotMat1stDerivatives [sigma ][l][dimA1];
45044629 }
45054630 }
45064631 }
@@ -4519,6 +4644,101 @@ void Mndo::RotateDiatomicTwoElecTwoCore2ndDerivativesToSpaceFrame(
45194644 // See Apendix in [DT_1977]
45204645 // Orbital mu and nu belong atom A,
45214646 // orbital lambda and sigma belong atomB.
4647+double Mndo::GetNddoRepulsionIntegralPointCharge(const Atom& atomA,
4648+ OrbitalType mu,
4649+ OrbitalType nu,
4650+ const Atom& atomB, // Point Charge
4651+ OrbitalType lambda,
4652+ OrbitalType sigma) const{
4653+ double value = 0.0;
4654+ double DA=0.0;
4655+ double DB=0.0;
4656+ double rhoA = 0.0;
4657+ double rhoB = 0.0;
4658+ double x = atomA.GetXyz()[0];
4659+ double y = atomA.GetXyz()[1];
4660+ double z = atomA.GetXyz()[2];
4661+ double rAB = sqrt(pow(x,2.0) + pow(y,2.0) + pow(z,2.0));
4662+ int lA = 0;
4663+ int lB = 0;
4664+ // (28) in [DT_1977]
4665+ if(mu == s && nu == s && lambda == s && sigma == s){
4666+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, sQ, sQ, rAB);
4667+ }
4668+ // (29) in [DT_1977]
4669+ else if(mu == s && nu == s && lambda == px && sigma == px){
4670+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, sQ, sQ, rAB);
4671+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, sQ, Qxx, rAB);
4672+ value = temp1 + temp2;
4673+ }
4674+ else if(mu == s && nu == s && lambda == py && sigma == py){
4675+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, sQ, sQ, rAB);
4676+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, sQ, Qyy, rAB);
4677+ value = temp1 + temp2;
4678+ }
4679+ // (30) in [DT_1977]
4680+ else if(mu == s && nu == s && lambda == pz && sigma == pz){
4681+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, sQ, sQ, rAB);
4682+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, sQ, Qzz, rAB);
4683+ value = temp1 + temp2;
4684+ }
4685+ // (31) in [DT_1977]
4686+ else if(mu == px && nu == px && lambda == s && sigma == s){
4687+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, sQ, sQ, rAB);
4688+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, Qxx, sQ, rAB);
4689+ value = temp1 + temp2;
4690+ }
4691+ else if(mu == py && nu == py && lambda == s && sigma == s){
4692+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, sQ, sQ, rAB);
4693+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, Qyy, sQ, rAB);
4694+ value = temp1 + temp2;
4695+ }
4696+ // (32) in [DT_1977]
4697+ else if(mu == pz && nu == pz && lambda == s && sigma == s){
4698+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, sQ, sQ, rAB);
4699+ double temp2 = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, Qzz, sQ, rAB);
4700+ value = temp1 + temp2;
4701+ }
4702+ // (38) in [DT_1977]
4703+ else if(mu == s && nu == pz && lambda == s && sigma == s){
4704+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, muz, sQ, rAB);
4705+ value = temp1;
4706+ }
4707+ else if(mu == pz && nu == s && lambda == s && sigma == s){
4708+ value = this->GetNddoRepulsionIntegralPointCharge(atomA, nu, mu, atomB, lambda, sigma);
4709+ }
4710+ // (41) in [DT_1977]
4711+ else if(mu == s && nu == s && lambda == s && sigma == pz){
4712+ double temp1 = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, sQ, muz, rAB);
4713+ value = temp1;
4714+ }
4715+ else if(mu == s && nu == s && lambda == pz && sigma == s){
4716+ value = this->GetNddoRepulsionIntegralPointCharge(atomA, mu, nu, atomB, sigma, lambda);
4717+ }
4718+ // d-orbitals
4719+ else if(mu == dxy || mu == dyz || mu == dzz || mu == dzx || mu == dxxyy ||
4720+ nu == dxy || nu == dyz || nu == dzz || nu == dzx || nu == dxxyy ||
4721+ lambda == dxy || lambda == dyz || lambda == dzz || lambda == dzx || lambda == dxxyy ||
4722+ sigma == dxy || sigma == dyz || sigma == dzz || sigma == dzx || sigma == dxxyy){
4723+
4724+ stringstream ss;
4725+ ss << this->errorMessageGetNddoRepulsionIntegral;
4726+ ss << this->errorMessageAtomA << AtomTypeStr(atomA.GetAtomType()) << endl;
4727+ ss << "\t" << this->errorMessageOrbitalType << OrbitalTypeStr(mu) << endl;
4728+ ss << "\t" << this->errorMessageOrbitalType << OrbitalTypeStr(nu) << endl;
4729+ ss << this->errorMessageAtomB << AtomTypeStr(atomB.GetAtomType()) << endl;
4730+ ss << "\t" << this->errorMessageOrbitalType << OrbitalTypeStr(lambda) << endl;
4731+ ss << "\t" << this->errorMessageOrbitalType << OrbitalTypeStr(sigma) << endl;
4732+ throw MolDSException(ss.str());
4733+ }
4734+ else{
4735+ value = 0.0;
4736+ }
4737+ return value;
4738+}
4739+// See Apendix in [DT_1977]
4740+// Orbital mu and nu belong atom A,
4741+// orbital lambda and sigma belong atomB.
45224742 double Mndo::GetNddoRepulsionIntegral(const Atom& atomA,
45234743 OrbitalType mu,
45244744 OrbitalType nu,
@@ -6228,6 +6448,230 @@ double Mndo::GetNddoRepulsionIntegral2ndDerivative(
62286448 }
62296449
62306450 // See Apendix in [DT_1977]
6451+double Mndo::GetSemiEmpiricalMultipoleInteractionPointCharge(const Atom& atomA,
6452+ const Atom& atomB,
6453+ MultipoleType multipoleA,
6454+ MultipoleType multipoleB,
6455+ double rAB) const{
6456+ double value = 0.0;
6457+ double DA = atomA.GetNddoDerivedParameterD(this->theory, multipoleA);
6458+ double DB = atomB.GetNddoDerivedParameterD(this->theory, multipoleB);
6459+ double rhoA = 0.0;
6460+ double rhoB = 0.0;
6461+ double a = rhoA + rhoB;
6462+
6463+ // Eq. (52) in [DT_1977]
6464+ if(multipoleA == sQ && multipoleB == sQ){
6465+ value = 1.0/sqrt(rAB*rAB + a*a);
6466+ }
6467+ // Eq. (53) in [DT_1977]
6468+ else if(multipoleA == sQ && multipoleB == muz){
6469+ double temp1 = ((rAB+DB)*(rAB+DB)) + (a*a);
6470+ double temp2 = ((rAB-DB)*(rAB-DB)) + (a*a);
6471+ value = 1.0/sqrt(temp1)/2.0 - 1.0/sqrt(temp2)/2.0;
6472+ }
6473+ else if(multipoleA == muz && multipoleB == sQ){
6474+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomB, atomA, multipoleB, multipoleA, rAB);
6475+ value *= -1.0;
6476+ }
6477+ // Eq. (54) in [DT_1977]
6478+ else if(multipoleA == sQ && multipoleB == Qxx){
6479+ double temp1 = (rAB*rAB) + (4.0*DB*DB) + (a*a);
6480+ double temp2 = (rAB*rAB) + (a*a);
6481+ value = 1.0/sqrt(temp1)/2.0 - 1.0/sqrt(temp2)/2.0;
6482+ }
6483+ else if(multipoleA == Qxx && multipoleB == sQ){
6484+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomB, atomA, multipoleB, multipoleA, rAB);
6485+ }
6486+ else if(multipoleA == sQ && multipoleB == Qyy){
6487+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, multipoleA, Qxx, rAB);
6488+ }
6489+ else if(multipoleA == Qyy && multipoleB == sQ){
6490+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomB, atomA, multipoleB, multipoleA, rAB);
6491+ }
6492+ // Eq. (55) in [DT_1977]
6493+ else if(multipoleA == sQ && multipoleB == Qzz){
6494+ double temp1 = ((rAB+2.0*DB)*(rAB+2.0*DB)) + (a*a);
6495+ double temp2 = (rAB*rAB) + (a*a);
6496+ double temp3 = ((rAB-2.0*DB)*(rAB-2.0*DB)) + (a*a);
6497+ value = 1.0/sqrt(temp1)/4.0 - 1.0/sqrt(temp2)/2.0 + 1.0/sqrt(temp3)/4.0;
6498+ }
6499+ else if(multipoleA == Qzz && multipoleB == sQ){
6500+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomB, atomA, multipoleB, multipoleA, rAB);
6501+ }
6502+ // Eq. (56) in [DT_1977]
6503+ else if(multipoleA == mux && multipoleB == mux){
6504+ double temp1 = (rAB*rAB) + ((DA-DB)*(DA-DB)) + (a*a);
6505+ double temp2 = (rAB*rAB) + ((DA+DB)*(DA+DB)) + (a*a);
6506+ value = 1.0/sqrt(temp1)/2.0 - 1.0/sqrt(temp2)/2.0;
6507+ }
6508+ else if(multipoleA == muy && multipoleB == muy){
6509+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, mux, mux, rAB);
6510+ }
6511+ // Eq. (57) in [DT_1977]
6512+ else if(multipoleA == muz && multipoleB == muz){
6513+ double temp1 = ((rAB+DA-DB)*(rAB+DA-DB)) + (a*a);
6514+ double temp2 = ((rAB+DA+DB)*(rAB+DA+DB)) + (a*a);
6515+ double temp3 = ((rAB-DA-DB)*(rAB-DA-DB)) + (a*a);
6516+ double temp4 = ((rAB-DA+DB)*(rAB-DA+DB)) + (a*a);
6517+ value = 1.0/sqrt(temp1)/4.0 - 1.0/sqrt(temp2)/4.0
6518+ -1.0/sqrt(temp3)/4.0 + 1.0/sqrt(temp4)/4.0;
6519+ }
6520+ // Eq. (58) in [DT_1977]
6521+ else if(multipoleA == mux && multipoleB == Qxz){
6522+ double temp1 = ((rAB-DB)*(rAB-DB)) + ((DA-DB)*(DA-DB)) + (a*a);
6523+ double temp2 = ((rAB-DB)*(rAB-DB)) + ((DA+DB)*(DA+DB)) + (a*a);
6524+ double temp3 = ((rAB+DB)*(rAB+DB)) + ((DA-DB)*(DA-DB)) + (a*a);
6525+ double temp4 = ((rAB+DB)*(rAB+DB)) + ((DA+DB)*(DA+DB)) + (a*a);
6526+ value =-1.0/sqrt(temp1)/4.0 + 1.0/sqrt(temp2)/4.0
6527+ +1.0/sqrt(temp3)/4.0 - 1.0/sqrt(temp4)/4.0;
6528+ }
6529+ else if(multipoleA == Qxz && multipoleB == mux){
6530+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomB, atomA, multipoleB, multipoleA, rAB);
6531+ value *= -1.0;
6532+ }
6533+ else if(multipoleA == muy && multipoleB == Qyz){
6534+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, mux, Qxz, rAB);
6535+ }
6536+ else if(multipoleA == Qyz && multipoleB == muy){
6537+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomB, atomA, multipoleB, multipoleA, rAB);
6538+ value *= -1.0;
6539+ }
6540+ // Eq. (59) in [DT_1977]
6541+ else if(multipoleA == muz && multipoleB == Qxx){
6542+ double temp1 = ((rAB+DA)*(rAB+DA)) + (4.0*DB*DB) + (a*a);
6543+ double temp2 = ((rAB-DA)*(rAB-DA)) + (4.0*DB*DB) + (a*a);
6544+ double temp3 = ((rAB+DA)*(rAB+DA)) + (a*a);
6545+ double temp4 = ((rAB-DA)*(rAB-DA)) + (a*a);
6546+ value =-1.0/sqrt(temp1)/4.0 + 1.0/sqrt(temp2)/4.0
6547+ +1.0/sqrt(temp3)/4.0 - 1.0/sqrt(temp4)/4.0;
6548+ }
6549+ else if(multipoleA == Qxx && multipoleB == muz){
6550+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomB, atomA, multipoleB, multipoleA, rAB);
6551+ value *= -1.0;
6552+ }
6553+ else if(multipoleA == muz && multipoleB == Qyy){
6554+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, muz, Qxx, rAB);
6555+ }
6556+ else if(multipoleA == Qyy && multipoleB == muz){
6557+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomB, atomA, multipoleB, multipoleA, rAB);
6558+ value *= -1.0;
6559+ }
6560+ // Eq. (60) in [DT_1977]
6561+ else if(multipoleA == muz && multipoleB == Qzz){
6562+ double temp1 = ((rAB+DA-2.0*DB)*(rAB+DA-2.0*DB)) + (a*a);
6563+ double temp2 = ((rAB-DA-2.0*DB)*(rAB-DA-2.0*DB)) + (a*a);
6564+ double temp3 = ((rAB+DA+2.0*DB)*(rAB+DA+2.0*DB)) + (a*a);
6565+ double temp4 = ((rAB-DA+2.0*DB)*(rAB-DA+2.0*DB)) + (a*a);
6566+ double temp5 = ((rAB+DA)*(rAB+DA)) + (a*a);
6567+ double temp6 = ((rAB-DA)*(rAB-DA)) + (a*a);
6568+ value =-1.0/sqrt(temp1)/8.0 + 1.0/sqrt(temp2)/8.0
6569+ -1.0/sqrt(temp3)/8.0 + 1.0/sqrt(temp4)/8.0
6570+ +1.0/sqrt(temp5)/4.0 - 1.0/sqrt(temp6)/4.0;
6571+ }
6572+ else if(multipoleA == Qzz && multipoleB == muz){
6573+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomB, atomA, multipoleB, multipoleA, rAB);
6574+ value *= -1.0;
6575+ }
6576+ // Eq. (61) in [DT_1977]
6577+ else if(multipoleA == Qxx && multipoleB == Qxx){
6578+ double temp1 = (rAB*rAB) + 4.0*((DA-DB)*(DA-DB)) + (a*a);
6579+ double temp2 = (rAB*rAB) + 4.0*((DA+DB)*(DA+DB)) + (a*a);
6580+ double temp3 = (rAB*rAB) + (4.0*DA*DA) + (a*a);
6581+ double temp4 = (rAB*rAB) + (4.0*DB*DB) + (a*a);
6582+ double temp5 = (rAB*rAB) + (a*a);
6583+ value = 1.0/sqrt(temp1)/8.0 + 1.0/sqrt(temp2)/8.0
6584+ -1.0/sqrt(temp3)/4.0 - 1.0/sqrt(temp4)/4.0
6585+ +1.0/sqrt(temp5)/4.0;
6586+ }
6587+ else if(multipoleA == Qyy && multipoleB == Qyy){
6588+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, Qxx, Qxx, rAB);
6589+ }
6590+ // Eq. (62) in [DT_1977]
6591+ else if(multipoleA == Qxx && multipoleB == Qyy){
6592+ double temp1 = (rAB*rAB) + (4.0*DA*DA) + (4.0*DB*DB)+ (a*a);
6593+ double temp2 = (rAB*rAB) + (4.0*DA*DA) + (a*a);
6594+ double temp3 = (rAB*rAB) + (4.0*DB*DB) + (a*a);
6595+ double temp4 = (rAB*rAB) + (a*a);
6596+ value = 1.0/sqrt(temp1)/4.0 - 1.0/sqrt(temp2)/4.0
6597+ -1.0/sqrt(temp3)/4.0 + 1.0/sqrt(temp4)/4.0;
6598+ }
6599+ else if(multipoleA == Qyy && multipoleB == Qxx){
6600+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomB, atomA, multipoleB, multipoleA, rAB);
6601+ }
6602+ // Eq. (63) in [DT_1977]
6603+ else if(multipoleA == Qxx && multipoleB == Qzz){
6604+ double temp1 = ((rAB-2.0*DB)*(rAB-2.0*DB)) + (4.0*DA*DA) + (a*a);
6605+ double temp2 = ((rAB+2.0*DB)*(rAB+2.0*DB)) + (4.0*DA*DA) + (a*a);
6606+ double temp3 = ((rAB-2.0*DB)*(rAB-2.0*DB)) + (a*a);
6607+ double temp4 = ((rAB+2.0*DB)*(rAB+2.0*DB)) + (a*a);
6608+ double temp5 = (rAB*rAB) + (4.0*DA*DA) + (a*a);
6609+ double temp6 = (rAB*rAB) + (a*a);
6610+ value = 1.0/sqrt(temp1)/8.0 + 1.0/sqrt(temp2)/8.0
6611+ -1.0/sqrt(temp3)/8.0 - 1.0/sqrt(temp4)/8.0
6612+ -1.0/sqrt(temp5)/4.0 + 1.0/sqrt(temp6)/4.0;
6613+ }
6614+ else if(multipoleA == Qzz && multipoleB == Qxx){
6615+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomB, atomA, multipoleB, multipoleA, rAB);
6616+ }
6617+ else if(multipoleA == Qyy && multipoleB == Qzz){
6618+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, Qxx, multipoleB, rAB);
6619+ }
6620+ else if(multipoleA == Qzz && multipoleB == Qyy){
6621+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomB, atomA, multipoleB, multipoleA, rAB);
6622+ }
6623+ // Eq. (64) in [DT_1977]
6624+ else if(multipoleA == Qzz && multipoleB == Qzz){
6625+ double temp1 = ((rAB+2.0*DA-2.0*DB)*(rAB+2.0*DA-2.0*DB)) + (a*a);
6626+ double temp2 = ((rAB+2.0*DA+2.0*DB)*(rAB+2.0*DA+2.0*DB)) + (a*a);
6627+ double temp3 = ((rAB-2.0*DA-2.0*DB)*(rAB-2.0*DA-2.0*DB)) + (a*a);
6628+ double temp4 = ((rAB-2.0*DA+2.0*DB)*(rAB-2.0*DA+2.0*DB)) + (a*a);
6629+ double temp5 = ((rAB+2.0*DA)*(rAB+2.0*DA)) + (a*a);
6630+ double temp6 = ((rAB-2.0*DA)*(rAB-2.0*DA)) + (a*a);
6631+ double temp7 = ((rAB+2.0*DB)*(rAB+2.0*DB)) + (a*a);
6632+ double temp8 = ((rAB-2.0*DB)*(rAB-2.0*DB)) + (a*a);
6633+ double temp9 = (rAB*rAB) + (a*a);
6634+ value = 1.0/sqrt(temp1)/16.0 + 1.0/sqrt(temp2)/16.0
6635+ +1.0/sqrt(temp3)/16.0 + 1.0/sqrt(temp4)/16.0
6636+ -1.0/sqrt(temp5)/8.0 - 1.0/sqrt(temp6)/8.0
6637+ -1.0/sqrt(temp7)/8.0 - 1.0/sqrt(temp8)/8.0
6638+ +1.0/sqrt(temp9)/4.0;
6639+ }
6640+ // Eq. (65) in [DT_1977]
6641+ else if(multipoleA == Qxz && multipoleB == Qxz){
6642+ double temp1 = ((rAB+DA-DB)*(rAB+DA-DB)) + ((DA-DB)*(DA-DB)) + (a*a);
6643+ double temp2 = ((rAB+DA-DB)*(rAB+DA-DB)) + ((DA+DB)*(DA+DB)) + (a*a);
6644+ double temp3 = ((rAB+DA+DB)*(rAB+DA+DB)) + ((DA-DB)*(DA-DB)) + (a*a);
6645+ double temp4 = ((rAB+DA+DB)*(rAB+DA+DB)) + ((DA+DB)*(DA+DB)) + (a*a);
6646+ double temp5 = ((rAB-DA-DB)*(rAB-DA-DB)) + ((DA-DB)*(DA-DB)) + (a*a);
6647+ double temp6 = ((rAB-DA-DB)*(rAB-DA-DB)) + ((DA+DB)*(DA+DB)) + (a*a);
6648+ double temp7 = ((rAB-DA+DB)*(rAB-DA+DB)) + ((DA-DB)*(DA-DB)) + (a*a);
6649+ double temp8 = ((rAB-DA+DB)*(rAB-DA+DB)) + ((DA+DB)*(DA+DB)) + (a*a);
6650+ value = 1.0/sqrt(temp1)/8.0 - 1.0/sqrt(temp2)/8.0
6651+ -1.0/sqrt(temp3)/8.0 + 1.0/sqrt(temp4)/8.0
6652+ -1.0/sqrt(temp5)/8.0 + 1.0/sqrt(temp6)/8.0
6653+ +1.0/sqrt(temp7)/8.0 - 1.0/sqrt(temp8)/8.0;
6654+ }
6655+ else if(multipoleA == Qyz && multipoleB == Qyz){
6656+ value = this->GetSemiEmpiricalMultipoleInteractionPointCharge(atomA, atomB, Qxz, Qxz, rAB);
6657+ }
6658+ // Eq. (66) in [DT_1977]
6659+ else if(multipoleA == Qxy && multipoleB == Qxy){
6660+ double temp1 = (rAB*rAB) + 2.0*((DA-DB)*(DA-DB)) + (a*a);
6661+ double temp2 = (rAB*rAB) + 2.0*((DA+DB)*(DA+DB)) + (a*a);
6662+ double temp3 = (rAB*rAB) + 2.0*(DA*DA) + 2.0*(DB*DB) + (a*a);
6663+ value = 1.0/sqrt(temp1)/4.0 + 1.0/sqrt(temp2)/4.0
6664+ -1.0/sqrt(temp3)/2.0;
6665+ }
6666+ else{
6667+ stringstream ss;
6668+ ss << this->errorMessageGetSemiEmpiricalMultipoleInteractionBadMultipoles;
6669+ ss << this->errorMessageMultipoleA << MultipoleTypeStr(multipoleA) << endl;
6670+ ss << this->errorMessageMultipoleB << MultipoleTypeStr(multipoleB) << endl;
6671+ throw MolDSException(ss.str());
6672+ }
6673+ return value;
6674+}
62316675 double Mndo::GetSemiEmpiricalMultipoleInteraction(const Atom& atomA,
62326676 const Atom& atomB,
62336677 MultipoleType multipoleA,
--- a/src/mndo/Mndo.h
+++ b/src/mndo/Mndo.h
@@ -36,13 +36,13 @@ protected:
3636 std::string errorMessageGetNddoRepulsionIntegral;
3737 std::string errorMessageGetNddoRepulsionIntegral1stDerivative;
3838 std::string errorMessageGetNddoRepulsionIntegral2ndDerivative;
39- std::string errorMessageCalcDiatomicTwoElecTwoCoreNullMatrix;
40- std::string errorMessageCalcTwoElecTwoCoreNullMatrix;
41- std::string errorMessageCalcDiatomicTwoElecTwoCoreSameAtoms;
42- std::string errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesSameAtoms;
43- std::string errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesSameAtoms;
44- std::string errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesNullMatrix;
45- std::string errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesNullMatrix;
39+ std::string errorMessageCalcDiatomicTwoElecsTwoCoresNullMatrix;
40+ std::string errorMessageCalcTwoElecsTwoCoresNullMatrix;
41+ std::string errorMessageCalcDiatomicTwoElecsTwoCoresSameAtoms;
42+ std::string errorMessageCalcDiatomicTwoElecsTwoCores1stDerivativesSameAtoms;
43+ std::string errorMessageCalcDiatomicTwoElecsTwoCores2ndDerivativesSameAtoms;
44+ std::string errorMessageCalcDiatomicTwoElecsTwoCores1stDerivativesNullMatrix;
45+ std::string errorMessageCalcDiatomicTwoElecsTwoCores2ndDerivativesNullMatrix;
4646 virtual void SetMessages();
4747 virtual void SetEnableAtomTypes();
4848 virtual void CalcSCFProperties();
@@ -63,7 +63,7 @@ protected:
6363 double const* const* gammaAB,
6464 double const* const* orbitalElectronPopulation,
6565 double const* atomicElectronPopulation,
66- double const* const* const* const* const* const* twoElecTwoCore,
66+ double const* const* const* const* const* const* twoElecsTwoAtomCores,
6767 bool isGuess) const;
6868 virtual double GetFockOffDiagElement(const MolDS_base_atoms::Atom& atomA,
6969 const MolDS_base_atoms::Atom& atomB,
@@ -74,7 +74,7 @@ protected:
7474 double const* const* gammaAB,
7575 double const* const* overelap,
7676 double const* const* orbitalElectronPopulation,
77- double const* const* const* const* const* const* twoElecTwoCore,
77+ double const* const* const* const* const* const* twoElecsTwoAtomCores,
7878 bool isGuess) const;
7979 virtual void CalcDiatomicOverlapAOsInDiatomicFrame(double** diatomicOverlapAOs,
8080 const MolDS_base_atoms::Atom& atomA,
@@ -91,8 +91,8 @@ protected:
9191 virtual double GetExchangeInt(MolDS_base::OrbitalType orbital1,
9292 MolDS_base::OrbitalType orbital2,
9393 const MolDS_base_atoms::Atom& atom) const;
94- virtual void CalcTwoElecTwoCore(double****** twoElecTwoCore,
95- const MolDS_base::Molecule& molecule) const;
94+ virtual void CalcTwoElecsTwoCores(double****** twoElecsTwoAtomCores,
95+ const MolDS_base::Molecule& molecule) const;
9696 virtual double GetMolecularIntegralElement(int moI,
9797 int moJ,
9898 int moK,
@@ -112,7 +112,7 @@ private:
112112 std::string errorMessageMultipoleB;
113113 std::string messageHeatsFormation;
114114 std::string messageHeatsFormationTitle;
115- double**** twoElecTwoCoreMpiBuff;
115+ double**** twoElecsTwoAtomCoresMpiBuff;
116116 double heatsFormation;
117117 double GetAuxiliaryDiatomCoreRepulsionEnergy(const MolDS_base_atoms::Atom& atomA,
118118 const MolDS_base_atoms::Atom& atomB,
@@ -136,8 +136,8 @@ private:
136136 double const* const* const* const* orbitalElectronPopulation1stDerivs,
137137 double const* const* const* const* diatomicOverlapAOs1stDerivs,
138138 double const* const* const* const* const* diatomicOverlapAOs2ndDerivs,
139- double const* const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs,
140- double const* const* const* const* const* const* const* diatomicTwoElecTwoCore2ndDerivs) const;
139+ double const* const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs,
140+ double const* const* const* const* const* const* const* diatomicTwoElecsTwoCores2ndDerivs) const;
141141 double GetHessianElementDifferentAtomsSCF(int indexAtomA,
142142 int indexAtomB,
143143 MolDS_base::CartesianType axisA,
@@ -146,18 +146,18 @@ private:
146146 double const* const* const* const* orbitalElectronPopulation1stDerivs,
147147 double const* const* const* const* diatomicOverlapAOs1stDerivs,
148148 double const* const* const* const* const* diatomicOverlapAOs2ndDerivs,
149- double const* const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs,
150- double const* const* const* const* const* const* const* diatomicTwoElecTwoCore2ndDerivs) const;
149+ double const* const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs,
150+ double const* const* const* const* const* const* const* diatomicTwoElecsTwoCores2ndDerivs) const;
151151 void MallocTempMatricesEachThreadCalcHessianSCF(double***** diatomicOverlapAOs1stDerivs,
152152 double****** diatomicOverlapAOs2ndDerivs,
153- double******* diatomicTwoElecTwoCore1stDerivs,
154- double******** diatomicTwoElecTwoCore2ndDerivs,
153+ double******* diatomicTwoElecsTwoCores1stDerivs,
154+ double******** diatomicTwoElecsTwoCores2ndDerivs,
155155 double*** tmpRotMat,
156156 double*** tmpRotMat1stDeriv,
157157 double**** tmpRotMat1stDerivs,
158158 double***** tmpRotMat2ndDerivs,
159- double***** tmpDiatomicTwoElecTwoCore,
160- double****** tmpDiatomicTwoElecTwoCore1stDerivs,
159+ double***** tmpDiatomicTwoElecsTwoCores,
160+ double****** tmpDiatomicTwoElecsTwoCores1stDerivs,
161161 double*** tmpDiaOverlapAOsInDiaFrame,
162162 double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
163163 double*** tmpDiaOverlapAOs2ndDerivInDiaFrame,
@@ -169,14 +169,14 @@ private:
169169 double** tmpVectorBC) const;
170170 void FreeTempMatricesEachThreadCalcHessianSCF(double***** diatomicOverlapAOs1stDerivs,
171171 double****** diatomicOverlapAOs2ndDerivs,
172- double******* diatomicTwoElecTwoCore1stDerivs,
173- double******** diatomicTwoElecTwoCore2ndDerivs,
172+ double******* diatomicTwoElecsTwoCores1stDerivs,
173+ double******** diatomicTwoElecsTwoCores2ndDerivs,
174174 double*** tmpRotMat,
175175 double*** tmpRotMat1stDeriv,
176176 double**** tmpRotMat1stDerivs,
177177 double***** tmpRotMat2ndDerivs,
178- double***** tmpDiatomicTwoElecTwoCore,
179- double****** tmpDiatomicTwoElecTwoCore1stDerivs,
178+ double***** tmpDiatomicTwoElecsTwoCores,
179+ double****** tmpDiatomicTwoElecsTwoCores1stDerivs,
180180 double*** tmpDiaOverlapAOsInDiaFrame,
181181 double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
182182 double*** tmpDiaOverlapAOs2ndDerivInDiaFrame,
@@ -193,7 +193,7 @@ private:
193193 MolDS_base::CartesianType axisA1,
194194 MolDS_base::CartesianType axisA2,
195195 double const* const* orbitalElectronPopulation,
196- double const* const* const* const* const* const* diatomicTwoElecTwoCore2ndDerivs) const;
196+ double const* const* const* const* const* const* diatomicTwoElecsTwoCores2ndDerivs) const;
197197 double GetAuxiliaryHessianElement2(int mu,
198198 int nu,
199199 int indexAtomA,
@@ -202,7 +202,7 @@ private:
202202 MolDS_base::CartesianType axisA,
203203 MolDS_base::CartesianType axisB,
204204 double const* const* const* const* orbitalElectronPopulation1stDerivs,
205- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs) const;
205+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const;
206206 double GetAuxiliaryHessianElement3(int lambda,
207207 int sigma,
208208 int indexAtomA,
@@ -210,7 +210,7 @@ private:
210210 MolDS_base::CartesianType axisA1,
211211 MolDS_base::CartesianType axisA2,
212212 double const* const* orbitalElectronPopulation,
213- double const* const* const* const* const* const* diatomicTwoElecTwoCore2ndDerivs) const;
213+ double const* const* const* const* const* const* diatomicTwoElecsTwoCores2ndDerivs) const;
214214 double GetAuxiliaryHessianElement4(int lambda,
215215 int sigma,
216216 int indexAtomA,
@@ -219,7 +219,7 @@ private:
219219 MolDS_base::CartesianType axisA,
220220 MolDS_base::CartesianType axisB,
221221 double const* const* const* const* orbitalElectronPopulation1stDerivs,
222- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs) const;
222+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const;
223223 double GetAuxiliaryHessianElement5(int mu,
224224 int lambda,
225225 int indexAtomA,
@@ -246,7 +246,7 @@ private:
246246 MolDS_base::CartesianType axisA1,
247247 MolDS_base::CartesianType axisA2,
248248 double const* const* orbitalElectronPopulation,
249- double const* const* const* const* const* const* diatomicTwoElecTwoCore2ndDerivs) const;
249+ double const* const* const* const* const* const* diatomicTwoElecsTwoCores2ndDerivs) const;
250250 double GetAuxiliaryHessianElement8(int mu,
251251 int nu,
252252 int lambda,
@@ -258,7 +258,7 @@ private:
258258 MolDS_base::CartesianType axisB,
259259 double const* const* orbitalElectronPopulation,
260260 double const* const* const* const* orbitalElectronPopulation1stDerivs,
261- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs) const;
261+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const;
262262 void CalcOrbitalElectronPopulation1stDerivatives(double**** orbitalElectronPopulation1stDerivatives) const;
263263 void SolveCPHF(double** solutionsCPHF,
264264 const std::vector<MoIndexPair>& nonRedundantQIndeces,
@@ -271,12 +271,12 @@ private:
271271 const std::vector<MoIndexPair>& redundantQIndeces,
272272 int indexAtomA,
273273 MolDS_base::CartesianType axisA) const;
274- void MallocTempMatricesStaticFirstOrderFock(double****** diatomicTwoElecTwoCore1stDeriv,
274+ void MallocTempMatricesStaticFirstOrderFock(double****** diatomicTwoElecsTwoCores1stDeriv,
275275 double**** diatomicOverlapAOs1stDeriv,
276276 double*** tmpRotMat,
277277 double**** tmpRotMat1stDerivs,
278278 double***** tmpDiatomicTwoElecTwo) const;
279- void FreeTempMatricesStaticFirstOrderFock(double****** diatomicTwoElecTwoCore1stDeriv,
279+ void FreeTempMatricesStaticFirstOrderFock(double****** diatomicTwoElecsTwoCores1stDeriv,
280280 double**** diatomicOverlapAOs1stDeriv,
281281 double*** tmpRotMat,
282282 double**** tmpRotMat1stDerivs,
@@ -294,50 +294,66 @@ private:
294294 int indexAtomB,
295295 int mu,
296296 int nu,
297- double const* const* const* const* const* const* twoElecTwoCore) const;
297+ double const* const* const* const* const* const* twoElecsTwoAtomCores) const;
298298 double GetElectronCoreAttraction1stDerivative(int indexAtomA,
299299 int indexAtomB,
300300 int mu,
301301 int nu,
302- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivatives,
302+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivatives,
303303 MolDS_base::CartesianType axisA) const;
304- void CalcDiatomicTwoElecTwoCore(double**** matrix,
304+ void CalcDiatomicTwoElecsTwoCoresPointCharge(double**** matrix,
305305 double* tmpVec,
306306 double** tmpRotMat,
307307 double** tmpMatrixBC,
308308 double* tmpVectorBC,
309- int indexAtomA,
310- int indexAtomB) const;
311- void CalcDiatomicTwoElecTwoCore1stDerivatives(double***** matrix,
312- double** tmpRotMat,
313- double*** tmpRotMat1stDerivs,
314- double**** tmpDiatomicTwoElecTwoCore,
315- int indexAtomA,
316- int indexAtomB) const;
317- void CalcDiatomicTwoElecTwoCore2ndDerivatives(double****** matrix,
318- double** tmpRotMat,
319- double*** tmpRotMat1stDerivs,
320- double**** tmpRotMat2ndDerivs,
321- double**** tmpDiatomicTwoElecTwoCore,
322- double***** tmpDiatomicTwoElecTwoCore1stDerivs,
323- int indexAtomA,
324- int indexAtomB) const;
325- void RotateDiatomicTwoElecTwoCoreToSpaceFrame(double**** matrix,
326- double* tmpVec,
327- double const* const* rotatingMatrix,
328- double** tmpMatrixBC,
329- double* tmpVectorBC) const;
330- void RotateDiatomicTwoElecTwoCore1stDerivativesToSpaceFrame(double***** matrix,
331- double const* const* const* const* diatomicTwoElecTwoCore,
332- double const* const* rotatingMatrix,
333- double const* const* const* rotMat1stDerivatives) const;
334- void RotateDiatomicTwoElecTwoCore2ndDerivativesToSpaceFrame(double****** matrix,
335- double const* const* const* const* diatomicTwoElecTwoCore,
336- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivatives,
337- double const* const* rotatingMatrix,
338- double const* const* const* rotMat1stDerivatives,
339- double const* const* const* const* rotMat2ndDerivatives) const;
340- void MallocTempMatricesRotateDiatomicTwoElecTwoCore1stDerivs(double*** twiceRotatingMatrix,
309+ const MolDS_base_atoms::Atom& atomA,
310+ const MolDS_base_atoms::Atom& pc) const;
311+ void CalcDiatomicTwoElecsTwoCores(double**** matrix,
312+ double* tmpVec,
313+ double** tmpRotMat,
314+ double** tmpMatrixBC,
315+ double* tmpVectorBC,
316+ int indexAtomA,
317+ int indexAtomB) const;
318+ void CalcDiatomicTwoElecsTwoCores1stDerivatives(double***** matrix,
319+ double** tmpRotMat,
320+ double*** tmpRotMat1stDerivs,
321+ double**** tmpDiatomicTwoElecsTwoCores,
322+ int indexAtomA,
323+ int indexAtomB) const;
324+ void CalcDiatomicTwoElecsTwoCores2ndDerivatives(double****** matrix,
325+ double** tmpRotMat,
326+ double*** tmpRotMat1stDerivs,
327+ double**** tmpRotMat2ndDerivs,
328+ double**** tmpDiatomicTwoElecsTwoCores,
329+ double***** tmpDiatomicTwoElecsTwoCores1stDerivs,
330+ int indexAtomA,
331+ int indexAtomB) const;
332+ void RotateDiatomicTwoElecsTwoCoresToSpaceFrame(double**** matrix,
333+ double* tmpVec,
334+ double const* const* rotatingMatrix,
335+ double** tmpMatrixBC,
336+ double* tmpVectorBC) const;
337+ void RotateDiatomicTwoElecsTwoCores1stDerivativesToSpaceFrame(double***** matrix,
338+ double const* const* const* const* diatomicTwoElecsTwoCores,
339+ double const* const* rotatingMatrix,
340+ double const* const* const* rotMat1stDerivatives) const;
341+ void RotateDiatomicTwoElecsTwoCores2ndDerivativesToSpaceFrame(double****** matrix,
342+ double const* const* const* const* diatomicTwoElecsTwoCores,
343+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivatives,
344+ double const* const* rotatingMatrix,
345+ double const* const* const* rotMat1stDerivatives,
346+ double const* const* const* const* rotMat2ndDerivatives) const;
347+ void MallocTempMatricesRotateDiatomicTwoElecsTwoCores1stDerivs(double*** twiceRotatingMatrix,
348+ double*** twiceRotatingMatrixDerivA,
349+ double*** twiceRotatingMatrixDerivB,
350+ double*** oldMatrix,
351+ double*** rotatedMatrix,
352+ double** tmpRotatedVec,
353+ double*** tmpMatrix,
354+ double** tmpVector,
355+ double*** ptrDiatomic) const;
356+ void FreeTempMatricesRotateDiatomicTwoElecsTwoCores1stDerivs(double*** twiceRotatingMatrix,
341357 double*** twiceRotatingMatrixDerivA,
342358 double*** twiceRotatingMatrixDerivB,
343359 double*** oldMatrix,
@@ -346,15 +362,12 @@ private:
346362 double*** tmpMatrix,
347363 double** tmpVector,
348364 double*** ptrDiatomic) const;
349- void FreeTempMatricesRotateDiatomicTwoElecTwoCore1stDerivs(double*** twiceRotatingMatrix,
350- double*** twiceRotatingMatrixDerivA,
351- double*** twiceRotatingMatrixDerivB,
352- double*** oldMatrix,
353- double*** rotatedMatrix,
354- double** tmpRotatedVec,
355- double*** tmpMatrix,
356- double** tmpVector,
357- double*** ptrDiatomic) const;
365+ double GetNddoRepulsionIntegralPointCharge(const MolDS_base_atoms::Atom& atomA,
366+ MolDS_base::OrbitalType mu,
367+ MolDS_base::OrbitalType nu,
368+ const MolDS_base_atoms::Atom& atomB,
369+ MolDS_base::OrbitalType lambda,
370+ MolDS_base::OrbitalType sigma) const;
358371 double GetNddoRepulsionIntegral(const MolDS_base_atoms::Atom& atomA,
359372 MolDS_base::OrbitalType mu,
360373 MolDS_base::OrbitalType nu,
@@ -376,6 +389,11 @@ private:
376389 MolDS_base::OrbitalType sigma,
377390 MolDS_base::CartesianType axisA1,
378391 MolDS_base::CartesianType axisA2) const;
392+ double GetSemiEmpiricalMultipoleInteractionPointCharge(const MolDS_base_atoms::Atom& atomA,
393+ const MolDS_base_atoms::Atom& atomB,
394+ MolDS_base::MultipoleType multipoleA,
395+ MolDS_base::MultipoleType multipoleB,
396+ double rAB) const;
379397 double GetSemiEmpiricalMultipoleInteraction(const MolDS_base_atoms::Atom& atomA,
380398 const MolDS_base_atoms::Atom& atomB,
381399 MolDS_base::MultipoleType multipoleA,
@@ -392,7 +410,7 @@ private:
392410 MolDS_base::MultipoleType multipoleB,
393411 double rAB) const;
394412 void MallocTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
395- double****** diatomicTwoElecTwoCore1stDerivs,
413+ double****** diatomicTwoElecsTwoCores1stDerivs,
396414 double*** tmpDiaOverlapAOsInDiaFrame,
397415 double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
398416 double*** tmpRotMat,
@@ -402,9 +420,9 @@ private:
402420 double** tmpRotatedDiatomicOverlapVec,
403421 double*** tmpMatrixBC,
404422 double** tmpVectorBC,
405- double***** tmpDiatomicTwoElecTwoCore) const;
423+ double***** tmpDiatomicTwoElecsTwoCores) const;
406424 void FreeTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
407- double****** diatomicTwoElecTwoCore1stDerivs,
425+ double****** diatomicTwoElecsTwoCores1stDerivs,
408426 double*** tmpDiaOverlapAOsInDiaFrame,
409427 double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
410428 double*** tmpRotMat,
@@ -414,11 +432,11 @@ private:
414432 double** tmpRotatedDiatomicOverlapVec,
415433 double*** tmpMatrixBC,
416434 double** tmpVectorBC,
417- double***** tmpDiatomicTwoElecTwoCore) const;
435+ double***** tmpDiatomicTwoElecsTwoCores) const;
418436 void CalcForceSCFElecCoreAttractionPart(double* force,
419437 int indexAtomA,
420438 int indexAtomB,
421- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs) const;
439+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const;
422440 void CalcForceSCFOverlapAOsPart(double* force,
423441 int indexAtomA,
424442 int indexAtomB,
@@ -426,22 +444,22 @@ private:
426444 void CalcForceSCFTwoElecPart(double* force,
427445 int indexAtomA,
428446 int indexAtomB,
429- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs) const;
447+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const;
430448 void CalcForceExcitedStaticPart(double* force,
431449 int elecStateIndex,
432450 int indexAtomA,
433451 int indexAtomB,
434- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs) const;
452+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const;
435453 void CalcForceExcitedElecCoreAttractionPart(double* force,
436454 int elecStateIndex,
437455 int indexAtomA,
438456 int indexAtomB,
439- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs) const;
457+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const;
440458 void CalcForceExcitedTwoElecPart(double* force,
441459 int elecStateIndex,
442460 int indexAtomA,
443461 int indexAtomB,
444- double const* const* const* const* const* diatomicTwoElecTwoCore1stDerivs) const;
462+ double const* const* const* const* const* diatomicTwoElecsTwoCores1stDerivs) const;
445463
446464 };
447465
--- a/src/pm3/Pm3.cpp
+++ b/src/pm3/Pm3.cpp
@@ -94,20 +94,20 @@ void Pm3::SetMessages(){
9494 = "Error in pm3::Pm3::GetNddoRepulsionIntegral1stDerivative: Bad orbital is set.\n";
9595 this->errorMessageGetNddoRepulsionIntegral2ndDerivative
9696 = "Error in pm3::Pm3::GetNddoRepulsionIntegral2ndDerivative: Bad orbital is set.\n";
97- this->errorMessageCalcTwoElecTwoCoreNullMatrix
98- = "Error in pm3::Pm3::CalcTwoElecTwoCore: The two elec two core matrix is NULL.\n";
99- this->errorMessageCalcDiatomicTwoElecTwoCoreSameAtoms
100- = "Error in pm3::Pm3::CalcDiatomicTwoElecTwoCore: Atom A and B is same.\n";
101- this->errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesSameAtoms
102- = "Error in pm3::Pm3::CalcDiatomicTwoElecTwoCore1stDerivatives: Atom A and B is same.\n";
103- this->errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesSameAtoms
104- = "Error in pm3::Pm3::CalcDiatomicTwoElecTwoCore2ndDerivatives: Atom A and B is same.\n";
105- this->errorMessageCalcDiatomicTwoElecTwoCoreNullMatrix
106- = "Error in pm3::Pm3::CalcDiatomicTwoElecTwoCore: The two elec two core diatomic matrix is NULL.\n";
107- this->errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesNullMatrix
108- = "Error in pm3::Pm3::CalcDiatomicTwoElecTwoCore1stDerivatives: The two elec two core diatomic matrix is NULL.\n";
109- this->errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesNullMatrix
110- = "Error in pm3::Pm3::CalcDiatomicTwoElecTwoCore2ndDerivatives: The two elec two core diatomic matrix is NULL.\n";
97+ this->errorMessageCalcTwoElecsTwoCoresNullMatrix
98+ = "Error in pm3::Pm3::CalcTwoElecsTwoCores: The two elec two core matrix is NULL.\n";
99+ this->errorMessageCalcDiatomicTwoElecsTwoCoresSameAtoms
100+ = "Error in pm3::Pm3::CalcDiatomicTwoElecsTwoCores: Atom A and B is same.\n";
101+ this->errorMessageCalcDiatomicTwoElecsTwoCores1stDerivativesSameAtoms
102+ = "Error in pm3::Pm3::CalcDiatomicTwoElecsTwoCores1stDerivatives: Atom A and B is same.\n";
103+ this->errorMessageCalcDiatomicTwoElecsTwoCores2ndDerivativesSameAtoms
104+ = "Error in pm3::Pm3::CalcDiatomicTwoElecsTwoCores2ndDerivatives: Atom A and B is same.\n";
105+ this->errorMessageCalcDiatomicTwoElecsTwoCoresNullMatrix
106+ = "Error in pm3::Pm3::CalcDiatomicTwoElecsTwoCores: The two elec two core diatomic matrix is NULL.\n";
107+ this->errorMessageCalcDiatomicTwoElecsTwoCores1stDerivativesNullMatrix
108+ = "Error in pm3::Pm3::CalcDiatomicTwoElecsTwoCores1stDerivatives: The two elec two core diatomic matrix is NULL.\n";
109+ this->errorMessageCalcDiatomicTwoElecsTwoCores2ndDerivativesNullMatrix
110+ = "Error in pm3::Pm3::CalcDiatomicTwoElecsTwoCores2ndDerivatives: The two elec two core diatomic matrix is NULL.\n";
111111 this->errorMessageGetElectronicEnergyEnergyNotCalculated
112112 = "Error in pm3::Pm3::GetElectronicEnergy: Set electronic state is not calculated by CIS.\n";
113113 this->errorMessageGetElectronicEnergyNULLCISEnergy
--- a/src/pm3/Pm3D.cpp
+++ b/src/pm3/Pm3D.cpp
@@ -95,20 +95,20 @@ void Pm3D::SetMessages(){
9595 = "Error in pm3::Pm3D::GetNddoRepulsionIntegral1stDerivative: Bad orbital is set.\n";
9696 this->errorMessageGetNddoRepulsionIntegral2ndDerivative
9797 = "Error in pm3::Pm3D::GetNddoRepulsionIntegral2ndDerivative: Bad orbital is set.\n";
98- this->errorMessageCalcTwoElecTwoCoreNullMatrix
99- = "Error in pm3::Pm3D::CalcTwoElecTwoCore: The two elec two core matrix is NULL.\n";
100- this->errorMessageCalcDiatomicTwoElecTwoCoreSameAtoms
101- = "Error in pm3::Pm3D::CalcDiatomicTwoElecTwoCore: Atom A and B is same.\n";
102- this->errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesSameAtoms
103- = "Error in pm3::Pm3D::CalcDiatomicTwoElecTwoCore1stDerivatives: Atom A and B is same.\n";
104- this->errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesSameAtoms
105- = "Error in pm3::Pm3D::CalcDiatomicTwoElecTwoCore2ndDerivatives: Atom A and B is same.\n";
106- this->errorMessageCalcDiatomicTwoElecTwoCoreNullMatrix
107- = "Error in pm3::Pm3D::CalcDiatomicTwoElecTwoCore: The two elec two core diatomic matrix is NULL.\n";
108- this->errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesNullMatrix
109- = "Error in pm3::Pm3D::CalcDiatomicTwoElecTwoCore1stDerivatives: The two elec two core diatomic matrix is NULL.\n";
110- this->errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesNullMatrix
111- = "Error in pm3::Pm3D::CalcDiatomicTwoElecTwoCore2ndDerivatives: The two elec two core diatomic matrix is NULL.\n";
98+ this->errorMessageCalcTwoElecsTwoCoresNullMatrix
99+ = "Error in pm3::Pm3D::CalcTwoElecsTwoCores: The two elec two core matrix is NULL.\n";
100+ this->errorMessageCalcDiatomicTwoElecsTwoCoresSameAtoms
101+ = "Error in pm3::Pm3D::CalcDiatomicTwoElecsTwoCores: Atom A and B is same.\n";
102+ this->errorMessageCalcDiatomicTwoElecsTwoCores1stDerivativesSameAtoms
103+ = "Error in pm3::Pm3D::CalcDiatomicTwoElecsTwoCores1stDerivatives: Atom A and B is same.\n";
104+ this->errorMessageCalcDiatomicTwoElecsTwoCores2ndDerivativesSameAtoms
105+ = "Error in pm3::Pm3D::CalcDiatomicTwoElecsTwoCores2ndDerivatives: Atom A and B is same.\n";
106+ this->errorMessageCalcDiatomicTwoElecsTwoCoresNullMatrix
107+ = "Error in pm3::Pm3D::CalcDiatomicTwoElecsTwoCores: The two elec two core diatomic matrix is NULL.\n";
108+ this->errorMessageCalcDiatomicTwoElecsTwoCores1stDerivativesNullMatrix
109+ = "Error in pm3::Pm3D::CalcDiatomicTwoElecsTwoCores1stDerivatives: The two elec two core diatomic matrix is NULL.\n";
110+ this->errorMessageCalcDiatomicTwoElecsTwoCores2ndDerivativesNullMatrix
111+ = "Error in pm3::Pm3D::CalcDiatomicTwoElecsTwoCores2ndDerivatives: The two elec two core diatomic matrix is NULL.\n";
112112 this->errorMessageGetElectronicEnergyEnergyNotCalculated
113113 = "Error in pm3::Pm3D::GetElectronicEnergy: Set electronic state is not calculated by CIS.\n";
114114 this->errorMessageGetElectronicEnergyNULLCISEnergy
--- a/src/pm3/Pm3Pddg.cpp
+++ b/src/pm3/Pm3Pddg.cpp
@@ -95,20 +95,20 @@ void Pm3Pddg::SetMessages(){
9595 = "Error in pm3::Pm3Pddg::GetNddoRepulsionIntegral1stDerivative: Bad orbital is set.\n";
9696 this->errorMessageGetNddoRepulsionIntegral2ndDerivative
9797 = "Error in pm3::Pm3Pddg::GetNddoRepulsionIntegral2ndDerivative: Bad orbital is set.\n";
98- this->errorMessageCalcTwoElecTwoCoreNullMatrix
99- = "Error in pm3::Pm3Pddg::CalcTwoElecTwoCore: The two elec two core matrix is NULL.\n";
100- this->errorMessageCalcDiatomicTwoElecTwoCoreSameAtoms
101- = "Error in pm3::Pm3Pddg::CalcDiatomicTwoElecTwoCore: Atom A and B is same.\n";
102- this->errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesSameAtoms
103- = "Error in pm3::Pm3Pddg::CalcDiatomicTwoElecTwoCore1stDerivatives: Atom A and B is same.\n";
104- this->errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesSameAtoms
105- = "Error in pm3::Pm3Pddg::CalcDiatomicTwoElecTwoCore2ndDerivatives: Atom A and B is same.\n";
106- this->errorMessageCalcDiatomicTwoElecTwoCoreNullMatrix
107- = "Error in pm3::Pm3Pddg::CalcDiatomicTwoElecTwoCore: The two elec two core diatomic matrix is NULL.\n";
108- this->errorMessageCalcDiatomicTwoElecTwoCore1stDerivativesNullMatrix
109- = "Error in pm3::Pm3Pddg::CalcDiatomicTwoElecTwoCore1stDerivatives: The two elec two core diatomic matrix is NULL.\n";
110- this->errorMessageCalcDiatomicTwoElecTwoCore2ndDerivativesNullMatrix
111- = "Error in pm3::Pm3Pddg::CalcDiatomicTwoElecTwoCore2ndDerivatives: The two elec two core diatomic matrix is NULL.\n";
98+ this->errorMessageCalcTwoElecsTwoCoresNullMatrix
99+ = "Error in pm3::Pm3Pddg::CalcTwoElecsTwoCores: The two elec two core matrix is NULL.\n";
100+ this->errorMessageCalcDiatomicTwoElecsTwoCoresSameAtoms
101+ = "Error in pm3::Pm3Pddg::CalcDiatomicTwoElecsTwoCores: Atom A and B is same.\n";
102+ this->errorMessageCalcDiatomicTwoElecsTwoCores1stDerivativesSameAtoms
103+ = "Error in pm3::Pm3Pddg::CalcDiatomicTwoElecsTwoCores1stDerivatives: Atom A and B is same.\n";
104+ this->errorMessageCalcDiatomicTwoElecsTwoCores2ndDerivativesSameAtoms
105+ = "Error in pm3::Pm3Pddg::CalcDiatomicTwoElecsTwoCores2ndDerivatives: Atom A and B is same.\n";
106+ this->errorMessageCalcDiatomicTwoElecsTwoCoresNullMatrix
107+ = "Error in pm3::Pm3Pddg::CalcDiatomicTwoElecsTwoCores: The two elec two core diatomic matrix is NULL.\n";
108+ this->errorMessageCalcDiatomicTwoElecsTwoCores1stDerivativesNullMatrix
109+ = "Error in pm3::Pm3Pddg::CalcDiatomicTwoElecsTwoCores1stDerivatives: The two elec two core diatomic matrix is NULL.\n";
110+ this->errorMessageCalcDiatomicTwoElecsTwoCores2ndDerivativesNullMatrix
111+ = "Error in pm3::Pm3Pddg::CalcDiatomicTwoElecsTwoCores2ndDerivatives: The two elec two core diatomic matrix is NULL.\n";
112112 this->errorMessageGetElectronicEnergyEnergyNotCalculated
113113 = "Error in pm3::Pm3Pddg::GetElectronicEnergy: Set electronic state is not calculated by CIS.\n";
114114 this->errorMessageGetElectronicEnergyNULLCISEnergy
--- a/src/zindo/ZindoS.cpp
+++ b/src/zindo/ZindoS.cpp
@@ -224,7 +224,7 @@ double ZindoS::GetFockDiagElement(const Atom& atomA,
224224 double const* const* gammaAB,
225225 double const* const* orbitalElectronPopulation,
226226 double const* atomicElectronPopulation,
227- double const* const* const* const* const* const* twoElecTwoCore,
227+ double const* const* const* const* const* const* twoElecsTwoAtomCores,
228228 bool isGuess) const{
229229 double value=0.0;
230230 int firstAOIndexA = atomA.GetFirstAOIndex();
@@ -289,7 +289,7 @@ double ZindoS::GetFockOffDiagElement(const Atom& atomA,
289289 double const* const* gammaAB,
290290 double const* const* overlapAOs,
291291 double const* const* orbitalElectronPopulation,
292- double const* const* const* const* const* const* twoElecTwoCore,
292+ double const* const* const* const* const* const* twoElecsTwoAtomCores,
293293 bool isGuess) const{
294294 double value = 0.0;
295295 OrbitalType orbitalMu = atomA.GetValence(mu-atomA.GetFirstAOIndex());
@@ -585,8 +585,8 @@ double ZindoS::GetExchangeInt(OrbitalType orbital1, OrbitalType orbital2, const
585585 return value;
586586 }
587587
588-void ZindoS::CalcTwoElecTwoCore(double****** twoElecTwoCore,
589- const Molecule& molecule) const{
588+void ZindoS::CalcTwoElecsTwoCores(double****** twoElecsTwoAtomCores,
589+ const Molecule& molecule) const{
590590 this->CalcNishimotoMatagaMatrix(this->nishimotoMatagaMatrix, molecule);
591591 }
592592
@@ -3620,9 +3620,9 @@ double ZindoS::GetAuxiliaryKNRKRElement(int moI, int moJ, int moK, int moL) cons
36203620 return value;
36213621 }
36223622
3623-void ZindoS::CalcDiatomicTwoElecTwoCore1stDerivatives(double*** matrix,
3624- int indexAtomA,
3625- int indexAtomB) const{
3623+void ZindoS::CalcDiatomicTwoElecsTwoCores1stDerivatives(double*** matrix,
3624+ int indexAtomA,
3625+ int indexAtomB) const{
36263626 const Atom& atomA = *molecule->GetAtom(indexAtomA);
36273627 const int firstAOIndexA = atomA.GetFirstAOIndex();
36283628 const int lastAOIndexA = atomA.GetLastAOIndex();
@@ -3665,7 +3665,7 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
36653665 #pragma omp parallel
36663666 {
36673667 double*** diatomicOverlapAOs1stDerivs = NULL;
3668- double*** diatomicTwoElecTwoCore1stDerivs = NULL;
3668+ double*** diatomicTwoElecsTwoCores1stDerivs = NULL;
36693669 double** tmpDiaOverlapAOsInDiaFrame = NULL; // diatomic overlapAOs in diatomic frame
36703670 double** tmpDiaOverlapAOs1stDerivInDiaFrame = NULL; // first derivative of the diaOverlapAOs. This derivative is related to the distance between two atoms.
36713671 double** tmpRotMat = NULL; // rotating Matrix from the diatomic frame to space fixed frame.
@@ -3677,7 +3677,7 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
36773677 double* tmpVectorBC = NULL; // used in dgemmm
36783678 try{
36793679 MallocTempMatricesCalcForce(&diatomicOverlapAOs1stDerivs,
3680- &diatomicTwoElecTwoCore1stDerivs,
3680+ &diatomicTwoElecsTwoCores1stDerivs,
36813681 &tmpDiaOverlapAOsInDiaFrame,
36823682 &tmpDiaOverlapAOs1stDerivInDiaFrame,
36833683 &tmpRotMat,
@@ -3710,7 +3710,7 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
37103710 atomB);
37113711
37123712 // calc. first derivative of two elec two core interaction by Nishimoto-Mataga
3713- this->CalcDiatomicTwoElecTwoCore1stDerivatives(diatomicTwoElecTwoCore1stDerivs, a, b);
3713+ this->CalcDiatomicTwoElecsTwoCores1stDerivatives(diatomicTwoElecsTwoCores1stDerivs, a, b);
37143714
37153715 double coreRepulsion [CartesianType_end] = {0.0,0.0,0.0};
37163716 double forceElecCoreAttPart[CartesianType_end] = {0.0,0.0,0.0};
@@ -3723,7 +3723,7 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
37233723 // electron core attraction part (ground state)
37243724 forceElecCoreAttPart[i] = ( atomA.GetCoreCharge()*atomicElectronPopulation[b]
37253725 +atomB.GetCoreCharge()*atomicElectronPopulation[a])
3726- *diatomicTwoElecTwoCore1stDerivs[s][s][i];
3726+ *diatomicTwoElecsTwoCores1stDerivs[s][s][i];
37273727 }
37283728 double forceOverlapAOsPart [CartesianType_end] = {0.0,0.0,0.0};
37293729 double forceTwoElecPart [CartesianType_end] = {0.0,0.0,0.0};
@@ -3742,7 +3742,7 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
37423742 forceTwoElecPart[i] += (this->orbitalElectronPopulation[mu][mu]
37433743 *this->orbitalElectronPopulation[nu][nu]
37443744 -0.5*pow(this->orbitalElectronPopulation[mu][nu],2.0))
3745- *diatomicTwoElecTwoCore1stDerivs[mu-firstAOIndexA][nu-firstAOIndexB][i];
3745+ *diatomicTwoElecsTwoCores1stDerivs[mu-firstAOIndexA][nu-firstAOIndexB][i];
37463746 }
37473747 }
37483748 }
@@ -3767,7 +3767,7 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
37673767 n,
37683768 a,
37693769 b,
3770- diatomicTwoElecTwoCore1stDerivs);
3770+ diatomicTwoElecsTwoCores1stDerivs);
37713771 // sum up contributions from static part (excited state)
37723772 #pragma omp critical
37733773 {
@@ -3785,7 +3785,7 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
37853785 n,
37863786 a,
37873787 b,
3788- diatomicTwoElecTwoCore1stDerivs);
3788+ diatomicTwoElecsTwoCores1stDerivs);
37893789 // overlapAOs part (excited states)
37903790 double forceExcitedOverlapAOsPart[CartesianType_end] = {0.0,0.0,0.0};
37913791 this->CalcForceExcitedOverlapAOsPart(forceExcitedOverlapAOsPart,
@@ -3799,7 +3799,7 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
37993799 n,
38003800 a,
38013801 b,
3802- diatomicTwoElecTwoCore1stDerivs);
3802+ diatomicTwoElecsTwoCores1stDerivs);
38033803 // sum up contributions from response part (excited state)
38043804 #pragma omp critical
38053805 {
@@ -3820,7 +3820,7 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
38203820 ex.Serialize(ompErrors);
38213821 }
38223822 FreeTempMatricesCalcForce(&diatomicOverlapAOs1stDerivs,
3823- &diatomicTwoElecTwoCore1stDerivs,
3823+ &diatomicTwoElecsTwoCores1stDerivs,
38243824 &tmpDiaOverlapAOsInDiaFrame,
38253825 &tmpDiaOverlapAOs1stDerivInDiaFrame,
38263826 &tmpRotMat,
@@ -3925,7 +3925,7 @@ void ZindoS::CalcForce(const vector<int>& elecStates){
39253925 }
39263926
39273927 void ZindoS::MallocTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
3928- double**** diatomicTwoElecTwoCore1stDerivs,
3928+ double**** diatomicTwoElecsTwoCores1stDerivs,
39293929 double*** tmpDiaOverlapAOsInDiaFrame,
39303930 double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
39313931 double*** tmpRotMat,
@@ -3936,7 +3936,7 @@ void ZindoS::MallocTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
39363936 double*** tmpMatrixBC,
39373937 double** tmpVectorBC) const{
39383938 MallocerFreer::GetInstance()->Malloc<double>(diatomicOverlapAOs1stDerivs, OrbitalType_end, OrbitalType_end, CartesianType_end);
3939- MallocerFreer::GetInstance()->Malloc<double>(diatomicTwoElecTwoCore1stDerivs, OrbitalType_end, OrbitalType_end, CartesianType_end);
3939+ MallocerFreer::GetInstance()->Malloc<double>(diatomicTwoElecsTwoCores1stDerivs, OrbitalType_end, OrbitalType_end, CartesianType_end);
39403940 MallocerFreer::GetInstance()->Malloc<double>(tmpDiaOverlapAOsInDiaFrame, OrbitalType_end, OrbitalType_end);
39413941 MallocerFreer::GetInstance()->Malloc<double>(tmpDiaOverlapAOs1stDerivInDiaFrame, OrbitalType_end, OrbitalType_end);
39423942 MallocerFreer::GetInstance()->Malloc<double>(tmpRotMat, OrbitalType_end, OrbitalType_end);
@@ -3949,7 +3949,7 @@ void ZindoS::MallocTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
39493949 }
39503950
39513951 void ZindoS::FreeTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
3952- double**** diatomicTwoElecTwoCore1stDerivs,
3952+ double**** diatomicTwoElecsTwoCores1stDerivs,
39533953 double*** tmpDiaOverlapAOsInDiaFrame,
39543954 double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
39553955 double*** tmpRotMat,
@@ -3960,7 +3960,7 @@ void ZindoS::FreeTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
39603960 double*** tmpMatrixBC,
39613961 double** tmpVectorBC) const{
39623962 MallocerFreer::GetInstance()->Free<double>(diatomicOverlapAOs1stDerivs, OrbitalType_end, OrbitalType_end, CartesianType_end);
3963- MallocerFreer::GetInstance()->Free<double>(diatomicTwoElecTwoCore1stDerivs, OrbitalType_end, OrbitalType_end, CartesianType_end);
3963+ MallocerFreer::GetInstance()->Free<double>(diatomicTwoElecsTwoCores1stDerivs, OrbitalType_end, OrbitalType_end, CartesianType_end);
39643964 MallocerFreer::GetInstance()->Free<double>(tmpDiaOverlapAOsInDiaFrame, OrbitalType_end, OrbitalType_end);
39653965 MallocerFreer::GetInstance()->Free<double>(tmpDiaOverlapAOs1stDerivInDiaFrame, OrbitalType_end, OrbitalType_end);
39663966 MallocerFreer::GetInstance()->Free<double>(tmpRotMat, OrbitalType_end, OrbitalType_end);
@@ -3976,7 +3976,7 @@ void ZindoS::CalcForceExcitedStaticPart(double* force,
39763976 int elecStateIndex,
39773977 int indexAtomA,
39783978 int indexAtomB,
3979- double const* const* const* diatomicTwoElecTwoCore1stDerivs) const{
3979+ double const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
39803980 const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
39813981 const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
39823982 int firstAOIndexA = atomA.GetFirstAOIndex();
@@ -3991,9 +3991,9 @@ void ZindoS::CalcForceExcitedStaticPart(double* force,
39913991 -1.0*this->etaMatrixForce[elecStateIndex][mu][lambda]
39923992 *this->etaMatrixForce[elecStateIndex][mu][lambda];
39933993 force[i] += temp
3994- *diatomicTwoElecTwoCore1stDerivs[mu-firstAOIndexA]
3995- [lambda-firstAOIndexB]
3996- [i];
3994+ *diatomicTwoElecsTwoCores1stDerivs[mu-firstAOIndexA]
3995+ [lambda-firstAOIndexB]
3996+ [i];
39973997 }
39983998 }
39993999 }
@@ -4003,7 +4003,7 @@ void ZindoS::CalcForceExcitedElecCoreAttractionPart(double* force,
40034003 int elecStateIndex,
40044004 int indexAtomA,
40054005 int indexAtomB,
4006- double const* const* const* diatomicTwoElecTwoCore1stDerivs) const{
4006+ double const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
40074007 const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
40084008 const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
40094009 int firstAOIndexA = atomA.GetFirstAOIndex();
@@ -4012,7 +4012,7 @@ void ZindoS::CalcForceExcitedElecCoreAttractionPart(double* force,
40124012 for(int i=0; i<CartesianType_end; i++){
40134013 force[i] += this->zMatrixForce[elecStateIndex][mu][mu]
40144014 *atomB.GetCoreCharge()
4015- *diatomicTwoElecTwoCore1stDerivs[mu-firstAOIndexA][s][i];
4015+ *diatomicTwoElecsTwoCores1stDerivs[mu-firstAOIndexA][s][i];
40164016 }
40174017 }
40184018 }
@@ -4051,7 +4051,7 @@ void ZindoS::CalcForceExcitedTwoElecPart(double* force,
40514051 int elecStateIndex,
40524052 int indexAtomA,
40534053 int indexAtomB,
4054- double const* const* const* diatomicTwoElecTwoCore1stDerivs) const{
4054+ double const* const* const* diatomicTwoElecsTwoCores1stDerivs) const{
40554055 const Atom& atomA = *this->molecule->GetAtom(indexAtomA);
40564056 const Atom& atomB = *this->molecule->GetAtom(indexAtomB);
40574057 int firstAOIndexA = atomA.GetFirstAOIndex();
@@ -4063,15 +4063,15 @@ void ZindoS::CalcForceExcitedTwoElecPart(double* force,
40634063 for(int i=0; i<CartesianType_end; i++){
40644064 force[i] -= this->zMatrixForce[elecStateIndex][mu][mu]
40654065 *this->orbitalElectronPopulation[lambda][lambda]
4066- *diatomicTwoElecTwoCore1stDerivs[mu-firstAOIndexA]
4067- [lambda-firstAOIndexB]
4068- [i];
4066+ *diatomicTwoElecsTwoCores1stDerivs[mu-firstAOIndexA]
4067+ [lambda-firstAOIndexB]
4068+ [i];
40694069 force[i] += 0.50
40704070 *this->zMatrixForce[elecStateIndex][mu][lambda]
40714071 *this->orbitalElectronPopulation[mu][lambda]
4072- *diatomicTwoElecTwoCore1stDerivs[mu-firstAOIndexA]
4073- [lambda-firstAOIndexB]
4074- [i];
4072+ *diatomicTwoElecsTwoCores1stDerivs[mu-firstAOIndexA]
4073+ [lambda-firstAOIndexB]
4074+ [i];
40754075 }
40764076 }
40774077 }
--- a/src/zindo/ZindoS.h
+++ b/src/zindo/ZindoS.h
@@ -70,7 +70,7 @@ protected:
7070 double const* const* gammaAB,
7171 double const* const* orbitalElectronPopulation,
7272 double const* atomicElectronPopulation,
73- double const* const* const* const* const* const* twoElecTwoCore,
73+ double const* const* const* const* const* const* twoElecsTwoAtomCores,
7474 bool isGuess) const;
7575 virtual double GetFockOffDiagElement(const MolDS_base_atoms::Atom& atomA,
7676 const MolDS_base_atoms::Atom& atomB,
@@ -81,7 +81,7 @@ protected:
8181 double const* const* gammaAB,
8282 double const* const* overelap,
8383 double const* const* orbitalElectronPopulation,
84- double const* const* const* const* const* const* twoElecTwoCore,
84+ double const* const* const* const* const* const* twoElecsTwoAtomCores,
8585 bool isGuess) const;
8686 virtual void CalcDiatomicOverlapAOsInDiatomicFrame(double** diatomicOverlapAOs,
8787 const MolDS_base_atoms::Atom& atomA,
@@ -98,8 +98,8 @@ protected:
9898 virtual double GetExchangeInt(MolDS_base::OrbitalType orbital1,
9999 MolDS_base::OrbitalType orbital2,
100100 const MolDS_base_atoms::Atom& atom) const; // Apendix in [BZ_1979]
101- virtual void CalcTwoElecTwoCore(double****** twoElecTwoCore,
102- const MolDS_base::Molecule& molecule) const;
101+ virtual void CalcTwoElecsTwoCores(double****** twoElecsTwoAtomCores,
102+ const MolDS_base::Molecule& molecule) const;
103103 virtual double GetMolecularIntegralElement(int moI,
104104 int moJ,
105105 int moK,
@@ -268,11 +268,11 @@ private:
268268 void FreeDavidsonRoopCISTemporaryMtrices(double*** interactionMatrix,
269269 int interactionMatrixDimension,
270270 double** interactionEigenEnergies) const;
271- void CalcDiatomicTwoElecTwoCore1stDerivatives(double*** matrix,
272- int indexAtomA,
273- int indexAtomB) const;
271+ void CalcDiatomicTwoElecsTwoCores1stDerivatives(double*** matrix,
272+ int indexAtomA,
273+ int indexAtomB) const;
274274 void MallocTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
275- double**** diatomicTwoElecTwoCore1stDerivs,
275+ double**** diatomicTwoElecsTwoCores1stDerivs,
276276 double*** tmpDiaOverlapAOsInDiaFrame,
277277 double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
278278 double*** tmpRotMat,
@@ -283,7 +283,7 @@ private:
283283 double*** tmpMatrixBC,
284284 double** tmpVectorBC) const;
285285 void FreeTempMatricesCalcForce(double**** diatomicOverlapAOs1stDerivs,
286- double**** diatomicTwoElecTwoCore1stDerivs,
286+ double**** diatomicTwoElecsTwoCores1stDerivs,
287287 double*** tmpDiaOverlapAOsInDiaFrame,
288288 double*** tmpDiaOverlapAOs1stDerivInDiaFrame,
289289 double*** tmpRotMat,
@@ -297,17 +297,17 @@ private:
297297 int elecStateIndex,
298298 int indexAtomA,
299299 int indexAtomB,
300- double const* const* const* diatomicTwoElecTwoCore1stDerivs) const;
300+ double const* const* const* diatomicTwoElecsTwoCores1stDerivs) const;
301301 void CalcForceExcitedElecCoreAttractionPart(double* force,
302302 int elecStateIndex,
303303 int indexAtomA,
304304 int indexAtomB,
305- double const* const* const* diatomicTwoElecTwoCore1stDerivs) const;
305+ double const* const* const* diatomicTwoElecsTwoCores1stDerivs) const;
306306 void CalcForceExcitedTwoElecPart(double* force,
307307 int elecStateIndex,
308308 int indexAtomA,
309309 int indexAtomB,
310- double const* const* const* diatomicTwoElecTwoCore1stDerivs) const;
310+ double const* const* const* diatomicTwoElecsTwoCores1stDerivs) const;
311311 void CheckZMatrixForce(const std::vector<int>& elecStates);
312312 void CheckEtaMatrixForce(const std::vector<int>& elecStates);
313313 double GetZMatrixForceElement(double const* y,