Commit: 1843479d8400c6f31dfdab20a3081c6ba9bb2fc5 - MolDS (git)

Commit MetaInfo

• Revision: 1843479d8400c6f31dfdab20a3081c6ba9bb2fc5 (tree)
• Time: 2013-12-04 11:54:29
• Author: Mikiya Fujii <mikiya.fujii@gmai...>
• Commiter: Mikiya Fujii

Log Message

Coulomb interaction between cores and EPCs is implemented. #32469

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

Change Summary

• modified: doc/README.txt
• modified: doc/refferences.txt
• modified: src/base/InputParser.cpp
• modified: src/base/InputParser.h
• modified: src/cndo/Cndo2.cpp
• modified: src/cndo/Cndo2.h
• modified: src/mndo/Mndo.cpp
• modified: src/mndo/Mndo.h

Incremental Difference

--- a/doc/README.txt
+++ b/doc/README.txt
@@ -339,11 +339,11 @@ HOW TO WRITE INPUT:
          MOPLOT_END
 
    <Environmental Point Charge(EPC) method>
-      Environmental point charge method is a simplified method of the QM/MM 
-      because the environmental point changes are treated as atoms in the MM region.
+      Environmental point charge method is a simplified method of the QM/MM,
+      namely the environmental point changes are treated as atoms in the MM region.
       The differences between the QM/MM and EPC are summarized below:
          - Electrostatic interaction between QM and MM region:
-            QM/MM: Electrostatic interaction is mutually added to QM and MM atoms.
+            QM/MM: Electrostatic interaction may be mutually added to QM and MM atoms.
             EPC  : Electrostatic field caused by the EPCs affects the QM region
                    although the each EPC is not affected by electrostatic field 
                    caused by the QM atoms and other EPCs. 
@@ -351,6 +351,8 @@ HOW TO WRITE INPUT:
          - Van der Waals interaction between QM and MM region:
             QM/MM: Included.
             EPC  : Not included.
+      In this EPC method, core-core replustion between QM and MM atoms is implemented with
+      the method I (simple coulomb interaction: qq/r) of ref [LRCL_2000].
       This EPC method can be used with MNDO-series (MNDO, AM1, AM1-D, PM3, PM3-D, and PDDG/PM3) only.
       To use this environmental point charges method, write EPC-directive. 
 
@@ -361,18 +363,19 @@ HOW TO WRITE INPUT:
       
       -options
        "the cartesian coordinates and charge" is only prepared.
-       Namely, each line should containe 4 doubles. 
+       Namely, each line should containe 4 doubles and a term. 
        The first three doubles are the cartesian coordinates of 
        each environmental point charge in angstrom unit.
-       The last double is the charge in atomic unit, 
+       The term is "charge". The last double following the term, "charge", 
+       is the charge in atomic unit, 
        e.g. -1 and 1 mean charge of an electron and a proton, respectively.
        Multiple setting of the environmental point charge is approvable, of course.
 
        E.g. 
          EPC 
-            0.0 0.0 0.0 -1
-            2.2 1.5 3.0 -1.5
-            0.0 2.0 5.0  0.3
+            0.0 0.0 0.0 charge -1   
+            2.2 1.5 3.0 charge -1.5
+            0.0 2.0 5.0 charge  0.3
          EPC_END
 
    <Frequencies (Normal modes analysis)>

--- a/doc/refferences.txt
+++ b/doc/refferences.txt
@@ -50,6 +50,7 @@
 [BFB_1997] M. A. Blanco, M. Fl{\'o}rez, and M. Bermejo, J. Mol. Struct. 419, 19 (1997)
 [BB_1998] E. Besalu, J. M. Bofill, Theor. Chem. Acc. 100, 265 (1998)
 [I_1998] IWANAMI-SHOTEN "IWANAMI RIKAGAKUJITEN", ISBN4-00-080090-6 (1988)
+[LRCL_2000] F. J. Luque, N. Reuter, A. Cartier, and M. F. R.-Lopez, J. Phys. Chem. A 104, 10923 (2000)
 [RCJ_2002] M. P. Repasky, J. Chandrasekhar, and W. L. Jorgensen, J. Comp. Chem. 23, 1601 (2002)
 [BGRJ_2003] I. T.-Brohman, C. R. W. Guimar\~{a}es, M. P. Repasky, and W. L. Jorgensen, J. Comp. Chem. 25, 138 (2003)
 [S_2004] J. J. P. Stewart, J. Mol. Model. 10, 155 (2004)

--- a/src/base/InputParser.cpp
+++ b/src/base/InputParser.cpp
@@ -251,8 +251,9 @@ void InputParser::SetMessages(){
    this->stringGeometryEnd = "geometry_end";
 
    // Environmental Point Charge
-   this->stringEpc    = "epc";
-   this->stringEpcEnd = "epc_end";
+   this->stringEpc       = "epc";
+   this->stringEpcEnd    = "epc_end";
+   this->stringEpcCharge = "charge";
 
    // SCF
    this->stringScf                 = "scf";
@@ -500,12 +501,16 @@ int InputParser::ParseEpcsConfiguration(Molecule* molecule, vector<string>* inpu
       double x      = atof((*inputTerms)[parseIndex+0].c_str()) * Parameters::GetInstance()->GetAngstrom2AU();
       double y      = atof((*inputTerms)[parseIndex+1].c_str()) * Parameters::GetInstance()->GetAngstrom2AU();
       double z      = atof((*inputTerms)[parseIndex+2].c_str()) * Parameters::GetInstance()->GetAngstrom2AU();
-      double charge = atof((*inputTerms)[parseIndex+3].c_str());
+      parseIndex += 3;
+      double charge = 0.0;
+      if((*inputTerms)[parseIndex].compare(this->stringEpcCharge) == 0){
+         charge = atof((*inputTerms)[parseIndex+1].c_str());
+         parseIndex += 2;
+      }
       AtomType atomType = EPC;
       int index = molecule->GetNumberEpcs();
       Atom* atom = AtomFactory::Create(atomType, index, x, y, z, charge);
       molecule->AddEpc(atom);
-      parseIndex += 4;
    }
    return parseIndex;
 }

--- a/src/base/InputParser.h
+++ b/src/base/InputParser.h
@@ -177,6 +177,7 @@ private:
    // EPC
    std::string stringEpc;
    std::string stringEpcEnd;
+   std::string stringEpcCharge;
    // SCF
    std::string stringScf;
    std::string stringScfEnd;

--- a/src/cndo/Cndo2.cpp
+++ b/src/cndo/Cndo2.cpp
@@ -71,9 +71,10 @@ Cndo2::Cndo2(){
    //protected variables
    this->molecule = NULL;
    this->theory = CNDO2;
-   this->coreRepulsionEnergy = 0.0;
-   this->vdWCorrectionEnergy = 0.0;
-   this->matrixCISdimension = 0;
+   this->coreRepulsionEnergy  = 0.0;
+   this->coreEpcCoulombEnergy = 0.0;
+   this->vdWCorrectionEnergy  = 0.0;
+   this->matrixCISdimension   = 0;
    this->fockMatrix = NULL;
    this->energiesMO = NULL;
    this->orbitalElectronPopulation    = NULL;
@@ -228,12 +229,16 @@ void Cndo2::SetMessages(){
    this->messageUnpairedAtoms      = "\tUnpaired electron population:";
    this->messageUnpairedAtomsTitle = "\t\t\t\t| k-th eigenstate | i-th atom | atom type | Unpaired electron population[a.u.]| \n";
    this->messageElecEnergy = "\tElectronic energy(SCF):";
-   this->messageNoteElecEnergy = "\tNote that this electronic energy includes core-repulsions.\n\n";
-   this->messageNoteElecEnergyVdW = "\tNote that this electronic energy includes core-repulsions and vdW correction.\n\n";
+   this->messageNoteElecEnergy       = "\tNote that this electronic energy includes core-repulsions.\n\n";
+   this->messageNoteElecEnergyVdW    = "\tNote that this electronic energy includes core-repulsions and vdW correction.\n\n";
+   this->messageNoteElecEnergyEpcVdW = "\tNote that this electronic energy includes core-repulsions, core-EPC coulomb, and vdW correction.\n\n";
+   this->messageNoteElecEnergyEpc    = "\tNote that this electronic energy includes core-repulsions and core-EPC coulomb.\n\n";
    this->messageElecEnergyTitle = "\t\t\t\t|   [a.u.]   |   [eV]   |\n";
    this->messageUnitSec = "[s].";
    this->messageCoreRepulsionTitle = "\t\t\t\t|   [a.u.]   |   [eV]   |\n";
    this->messageCoreRepulsion = "\tCore repulsion energy:";
+   this->messageCoreEpcCoulombTitle = "\t\t\t\t\t\t\t\t|   [a.u.]   |   [eV]   |\n";
+   this->messageCoreEpcCoulomb = "\tCoulomb interaction between cores and EPCs energy:";
    this->messageVdWCorrectionTitle = "\t\t\t\t\t\t|   [a.u.]   |   [eV]   |\n";
    this->messageVdWCorrection = "\tEmpirical van der Waals correction:";
    this->messageElectronicDipoleMomentTitle = "\t\t\t\t\t|  x[a.u.]  |  y[a.u.]  |  z[a.u.]  |  magnitude[a.u.]  |\t\t|  x[debye]  |  y[debye]  |  z[debye]  |  magnitude[debye]  |\n";
@@ -341,6 +346,7 @@ void Cndo2::CheckEnableAtomType(const Molecule& molecule) const{
 }
 
 void Cndo2::CalcCoreRepulsionEnergy(){
+   // interaction between atoms
    double energy = 0.0;
    for(int i=0; i<this->molecule->GetNumberAtoms(); i++){
       for(int j=i+1; j<this->molecule->GetNumberAtoms(); j++){
@@ -348,6 +354,17 @@ void Cndo2::CalcCoreRepulsionEnergy(){
       }
    }
    this->coreRepulsionEnergy = energy;
+
+   // interaction between atoms and epcs
+   if(this->molecule->GetNumberEpcs()<=0){return;}
+   energy = 0.0;
+   for(int i=0; i<this->molecule->GetNumberAtoms(); i++){
+      for(int j=0; j<this->molecule->GetNumberEpcs(); j++){
+         energy += this->GetAtomCoreEpcCoulombEnergy(i, j);
+      }
+   }
+   this->coreEpcCoulombEnergy = energy;
+
 }
 
 double Cndo2::GetDiatomCoreRepulsionEnergy(int indexAtomA, int indexAtomB) const{
@@ -357,6 +374,11 @@ double Cndo2::GetDiatomCoreRepulsionEnergy(int indexAtomA, int indexAtomB) const
    return atomA.GetCoreCharge()*atomB.GetCoreCharge()/distance; 
 }
 
+double Cndo2::GetAtomCoreEpcCoulombEnergy(int indexAtom, int indexEpc) const{
+   // do nothiing
+   return 0.0;
+}
+
 // First derivative of diatomic core repulsion energy.
 // This derivative is related to the coordinate of atomA.
 double Cndo2::GetDiatomCoreRepulsion1stDerivative(int indexAtomA, int indexAtomB, 
@@ -652,12 +674,13 @@ void Cndo2::CalcSCFProperties(){
       this->CalcVdWCorrectionEnergy();
    }
    this->CalcElecSCFEnergy(&this->elecSCFEnergy, 
-                          *this->molecule, 
-                          this->energiesMO, 
-                          this->fockMatrix, 
-                          this->gammaAB,
-                          this->coreRepulsionEnergy,
-                          this->vdWCorrectionEnergy);
+                           *this->molecule, 
+                           this->energiesMO, 
+                           this->fockMatrix, 
+                           this->gammaAB,
+                           this->coreRepulsionEnergy,
+                           this->coreEpcCoulombEnergy,
+                           this->vdWCorrectionEnergy);
    this->CalcCoreDipoleMoment(this->coreDipoleMoment, *this->molecule);
    this->CalcElectronicDipoleMomentGroundState(this->electronicTransitionDipoleMoments, 
                                                this->cartesianMatrix,
@@ -990,9 +1013,15 @@ void Cndo2::OutputSCFEnergies() const{
    this->OutputLog(boost::format("%s\t%e\t%e\n") % this->messageElecEnergy
                                                  % this->elecSCFEnergy
                                                  % (this->elecSCFEnergy/eV2AU));
-   if(Parameters::GetInstance()->RequiresVdWSCF()){
+   if(Parameters::GetInstance()->RequiresVdWSCF() && this->molecule->GetNumberEpcs()<=0){
       this->OutputLog(this->messageNoteElecEnergyVdW);
    }
+   else if(Parameters::GetInstance()->RequiresVdWSCF() && 0<this->molecule->GetNumberEpcs()){
+      this->OutputLog(this->messageNoteElecEnergyEpcVdW);
+   }
+   else if(!Parameters::GetInstance()->RequiresVdWSCF() && 0<this->molecule->GetNumberEpcs()){
+      this->OutputLog(this->messageNoteElecEnergyEpc);
+   }
    else{
       this->OutputLog(this->messageNoteElecEnergy);
    }
@@ -1003,6 +1032,14 @@ void Cndo2::OutputSCFEnergies() const{
                                                    % this->coreRepulsionEnergy 
                                                    % (this->coreRepulsionEnergy/eV2AU));
 
+   // output coulomb interaction between atoms and epcs
+   if(0<this->molecule->GetNumberEpcs()){
+      this->OutputLog(this->messageCoreEpcCoulombTitle);
+      this->OutputLog(boost::format("%s\t%e\t%e\n\n") % this->messageCoreEpcCoulomb
+                                                      % this->coreEpcCoulombEnergy 
+                                                      % (this->coreEpcCoulombEnergy/eV2AU));
+   }
+
    // output van der Waals correction 
    if(Parameters::GetInstance()->RequiresVdWSCF()){
       this->OutputLog(this->messageVdWCorrectionTitle);
@@ -1194,6 +1231,7 @@ void Cndo2::CalcElecSCFEnergy(double* elecSCFEnergy,
                              double const* const* fockMatrix, 
                              double const* const* gammaAB, 
                              double coreRepulsionEnergy,
+                             double coreEpcCoulombEnergy,
                              double vdWCorrectionEnergy) const{
    double electronicEnergy = 0.0;
    // use density matrix for electronic energy
@@ -1273,7 +1311,10 @@ void Cndo2::CalcElecSCFEnergy(double* elecSCFEnergy,
    }
    */
 
-   *elecSCFEnergy = electronicEnergy + coreRepulsionEnergy + vdWCorrectionEnergy;
+   *elecSCFEnergy = electronicEnergy
+                   +coreRepulsionEnergy
+                   +coreEpcCoulombEnergy
+                   +vdWCorrectionEnergy;
 }
 
 void Cndo2::FreeElecEnergyMatrices(double*** fMatrix, 

--- a/src/cndo/Cndo2.h
+++ b/src/cndo/Cndo2.h
@@ -94,6 +94,7 @@ protected:
    MolDS_base::Molecule* molecule;
    MolDS_base::TheoryType theory;
    double       coreRepulsionEnergy;
+   double       coreEpcCoulombEnergy;
    double       vdWCorrectionEnergy;
    int          matrixCISdimension;
    double**     fockMatrix;
@@ -132,6 +133,7 @@ protected:
    double GetBondingAdjustParameterK(MolDS_base::ShellType shellA, 
                                      MolDS_base::ShellType shellB) const;
    virtual double GetDiatomCoreRepulsionEnergy(int indexAtomA, int indexAtomB) const;
+   virtual double GetAtomCoreEpcCoulombEnergy (int indexAtom,  int indexEpc) const;
    virtual double GetDiatomCoreRepulsion1stDerivative(int indexAtomA, 
                                                       int indexAtomB, 
                                                       MolDS_base::CartesianType axisA) const;
@@ -305,11 +307,15 @@ private:
    std::string messageElecEnergy;
    std::string messageNoteElecEnergy;
    std::string messageNoteElecEnergyVdW;
+   std::string messageNoteElecEnergyEpcVdW;
+   std::string messageNoteElecEnergyEpc;
    std::string messageElecEnergyTitle;
    std::string messageOcc;
    std::string messageUnOcc;
    std::string messageCoreRepulsionTitle;
    std::string messageCoreRepulsion;
+   std::string messageCoreEpcCoulombTitle;
+   std::string messageCoreEpcCoulomb;
    std::string messageVdWCorrectionTitle;
    std::string messageVdWCorrection;
    std::string messageElectronicDipoleMomentTitle;
@@ -512,6 +518,7 @@ private:
                           double const* const* fockMatrix, 
                           double const* const* gammaAB, 
                           double coreRepulsionEnergy,
+                          double coreEpcCoulombEnergy,
                           double vdWCorrectionEnergy) const;
    void FreeElecEnergyMatrices(double*** fMatrix, 
                                double*** hMatrix, 

--- a/src/mndo/Mndo.cpp
+++ b/src/mndo/Mndo.cpp
@@ -326,6 +326,13 @@ double Mndo::GetDiatomCoreRepulsionEnergy(int indexAtomA, int indexAtomB) const{
          *tmp;
 }
 
+double Mndo::GetAtomCoreEpcCoulombEnergy(int indexAtom, int indexEpc) const{
+   const Atom& atom = *this->molecule->GetAtom(indexAtom);
+   const Atom& epc  = *this->molecule->GetAtom(indexEpc);
+   double distance = this->molecule->GetDistanceAtomEpc(indexAtom, indexEpc);
+   return atom.GetCoreCharge()*epc.GetCoreCharge()/distance; 
+}
+
 // First derivative of diatomic core repulsion energy.
 // This derivative is related to the coordinate of atomA.
 double Mndo::GetDiatomCoreRepulsion1stDerivative(int indexAtomA,

--- a/src/mndo/Mndo.h
+++ b/src/mndo/Mndo.h
@@ -53,6 +53,7 @@ protected:
    virtual void CalcNormalModes(double** normalModes, double* normalForceConstants, const MolDS_base::Molecule& molecule) const;
    virtual void CalcForce(const std::vector<int>& elecStates);
    virtual double GetDiatomCoreRepulsionEnergy(int indexAtomA, int indexAtomB) const;
+   virtual double GetAtomCoreEpcCoulombEnergy (int indexAtom,  int indexEpc  ) const;
    virtual double GetDiatomCoreRepulsion1stDerivative(int indexAtomA,
                                                       int indexAtomB, 
                                                       MolDS_base::CartesianType axisA) const;