提交: 23e25438a28a981d14c37ce2fd98c04e36378031 - MolDS (git)

Commit MetaInfo

• 修訂: 23e25438a28a981d14c37ce2fd98c04e36378031 (tree)
• 時間: 2014-04-29 14:02:11
• 作者: Mikiya Fujii <mikiya.fujii@gmai...>
• Commiter: Mikiya Fujii

Log Message

Check for adding bare VdW is added. #33724

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

Change Summary

• modified: doc/README.txt
• modified: src/am1/Am1.cpp
• modified: src/am1/Am1D.cpp
• modified: src/cndo/Cndo2.cpp
• modified: src/cndo/Cndo2.h
• modified: src/indo/Indo.cpp
• modified: src/mndo/Mndo.cpp
• modified: src/pm3/Pm3.cpp
• modified: src/pm3/Pm3D.cpp
• modified: src/pm3/Pm3Pddg.cpp
• modified: src/zindo/ZindoS.cpp

差異

--- a/doc/README.txt
+++ b/doc/README.txt
@@ -245,9 +245,10 @@ HOW TO WRITE INPUT:
        The default value of the "diis_end_error" is 10**(-8.0).
 
        "vdW" should be set as "yes" or "no". 
-       When "yes" is set, Grimmes's empirical van der Waals correction([G_2004]) is applied.
+       When "yes" is set, Grimmes's empirical van der Waals correction(D1, [G_2004]) is applied.
        Note that this empirical van der Waals correction is applied to the semiempirical theories 
        of which semiempirical parameters are not modified.
+       This "vdw" option can be used for H, C, N, O, F, S, and Cl.
        If user wants to use PM3-D or AM1-D of which semiempirical parameters are modified to be suite for vdW, 
        set theory-directive as "PM3-D" or "AM1-D". 
        When PM3-D or AM1-D is used, users do not need to set "vdW", "vdW_s6", and "vdW_d".
@@ -256,6 +257,14 @@ HOW TO WRITE INPUT:
        The default value of the "vdW" with the theories except for PM3-D and AM1-D is "no". 
        For PM3-D and AM1-D, this "vdW" is always "yes" whethere user sets or not.
 
+       "vdW_s6" is a scaling factor in the Grimme's van der Waals correction([G_2004]).
+       The default value of the "vdW_s6" is 1.4. 
+       For PM3-D and AM1-D, this "vdW_s6" is forced to be set as 1.4.
+
+       "vdW_d" is a damping factor in the Grimme's van der Waals correction([G_2004]).
+       The default value of the "vdW_d" is 23.0.
+       For PM3-D and AM1-D, this "vdW_s6" is forced to be set as 23.0.
+
        "sum_charges" is an option to calculate of summation of atomic charges in the ground state.
        To use this option, write
        "sum_charges first_atom_index last_atom_index"
@@ -265,14 +274,6 @@ HOW TO WRITE INPUT:
        Multiple setting of this "sum_charges" option is approvable, of course.
        If you want to calculate summation, same "sum_charges" option is available in CIS-directive.
 
-       "vdW_s6" is a scaling factor in the Grimme's van der Waals correction([G_2004]).
-       The default value of the "vdW_s6" is 1.4. 
-       For PM3-D and AM1-D, this "vdW_s6" is forced to be set as 1.4.
-
-       "vdW_d" is a damping factor in the Grimme's van der Waals correction([G_2004]).
-       The default value of the "vdW_d" is 23.0.
-       For PM3-D and AM1-D, this "vdW_s6" is forced to be set as 23.0.
-
        E.g.
          SCF
             max_iter 200

--- a/src/am1/Am1.cpp
+++ b/src/am1/Am1.cpp
@@ -77,6 +77,8 @@ void Am1::SetMessages(){
       = "Error in am1::Am1::SetMolecule: Total number of valence electrons is odd. totalNumberValenceElectrons=";
    this->errorMessageNotEnebleAtomType  
       = "Error in am1::Am1::CheckEnableAtomType: Non available atom is contained.\n";
+   this->errorMessageNotEnebleAtomTypeVdW
+      = "Error in am1::Am1::CheckEnableAtomTypeVdW: Non available atom to add VdW correction is contained.\n";
    this->errorMessageCoulombInt = "Error in base_am1::Am1::GetCoulombInt: Invalid orbitalType.\n";
    this->errorMessageExchangeInt = "Error in base_am1::Am1::GetExchangeInt: Invalid orbitalType.\n";
    this->errorMessageCalcCISMatrix

--- a/src/am1/Am1D.cpp
+++ b/src/am1/Am1D.cpp
@@ -79,6 +79,8 @@ void Am1D::SetMessages(){
       = "Error in am1::Am1D::SetMolecule: Total number of valence electrons is odd. totalNumberValenceElectrons=";
    this->errorMessageNotEnebleAtomType  
       = "Error in am1::Am1D::CheckEnableAtomType: Non available atom is contained.\n";
+   this->errorMessageNotEnebleAtomTypeVdW
+      = "Error in am1::Am1D::CheckEnableAtomTypeVdW: Non available atom to add VdW correction is contained.\n";
    this->errorMessageCoulombInt = "Error in am1::Am1D::GetCoulombInt: Invalid orbitalType.\n";
    this->errorMessageExchangeInt = "Error in am1::Am1D::GetExchangeInt: Invalid orbitalType.\n";
    this->errorMessageCalcCISMatrix

--- a/src/cndo/Cndo2.cpp
+++ b/src/cndo/Cndo2.cpp
@@ -101,6 +101,7 @@ Cndo2::Cndo2(){
    //protected methods
    this->SetMessages();
    this->SetEnableAtomTypes();
+   this->SetEnableAtomTypesVdW();
 
    //private variables
    this->elecSCFEnergy = 0.0;
@@ -152,7 +153,9 @@ void Cndo2::SetMessages(){
    this->errorMessageOddTotalValenceElectrions 
       = "Error in cndo::Cndo2::SetMolecule: Total number of valence electrons is odd. totalNumberValenceElectrons=";
    this->errorMessageNotEnebleAtomType  
-      = "Error in cndo::Cndo2::ChecEnableAtomType: Non available atom is contained.\n";
+      = "Error in cndo::Cndo2::CheckEnableAtomType: Non available atom is contained.\n";
+   this->errorMessageNotEnebleAtomTypeVdW
+      = "Error in cndo::Cndo2::CheckEnableAtomTypeVdW: Non available atom to add VdW correction is contained.\n";
    this->errorMessageAtomA = "Atom A is:\n";
    this->errorMessageAtomB = "Atom B is:\n";
    this->errorMessageAtomType = "\tatom type = ";
@@ -281,6 +284,17 @@ void Cndo2::SetEnableAtomTypes(){
    this->enableAtomTypes.push_back(Cl);
 }
 
+void Cndo2::SetEnableAtomTypesVdW(){
+   this->enableAtomTypesVdW.clear();
+   this->enableAtomTypesVdW.push_back(H);
+   this->enableAtomTypesVdW.push_back(C);
+   this->enableAtomTypesVdW.push_back(N);
+   this->enableAtomTypesVdW.push_back(O);
+   this->enableAtomTypesVdW.push_back(F);
+   this->enableAtomTypesVdW.push_back(S);
+   this->enableAtomTypesVdW.push_back(Cl);
+}
+
 TheoryType Cndo2::GetTheoryType() const{
    return this->theory;
 }
@@ -289,6 +303,7 @@ void Cndo2::SetMolecule(Molecule* molecule){
    this->molecule = molecule;
    this->CheckNumberValenceElectrons(*molecule);
    this->CheckEnableAtomType(*molecule);
+   this->CheckEnableAtomTypeVdW(*molecule);
 
    // malloc
    MallocerFreer::GetInstance()->Malloc<double>(&this->fockMatrix,
@@ -352,6 +367,34 @@ void Cndo2::CheckEnableAtomType(const Molecule& molecule) const{
    }
 }
 
+void Cndo2::CheckEnableAtomTypeVdW(const Molecule& molecule) const{
+   if(Parameters::GetInstance()->RequiresVdWSCF()){
+      if(this->theory == AM1D || this->theory == PM3D){
+         return;
+      }
+   }
+   else{
+      return;
+   }
+   
+   for(int i=0; i<molecule.GetAtomVect().size(); i++){
+      AtomType atomType = molecule.GetAtomVect()[i]->GetAtomType();
+      bool enable = false;
+      for(int j=0; j<this->enableAtomTypesVdW.size(); j++){
+         if(atomType == this->enableAtomTypesVdW[j]){
+            enable = true;
+            break;
+         }
+      }
+      if(!enable){
+         stringstream ss;
+         ss << this->errorMessageNotEnebleAtomTypeVdW;
+         ss << this->errorMessageAtomType << AtomTypeStr(atomType) << endl;
+         throw MolDSException(ss.str());
+      }
+   }
+}
+
 void Cndo2::CalcCoreRepulsionEnergy(){
    // interaction between atoms
    double energy = 0.0;
@@ -423,13 +466,13 @@ void Cndo2::CalcVdWCorrectionEnergy(){
 }
 
 // See damping function in (2) in [G_2004] ((11) in [G_2006])
-double Cndo2::GetVdwDampingValue(double vdWDistance, double distance) const{
+double Cndo2::GetVdWDampingValue(double vdWDistance, double distance) const{
    double dampingFactor = Parameters::GetInstance()->GetVdWDampingFactorSCF();
    return 1.0/(1.0+exp(-1.0*dampingFactor*(distance/vdWDistance - 1.0)));
 }
 
 // See damping function in (2) in [G_2004] ((11) in [G_2006])
-double Cndo2::GetVdwDampingValue1stDerivative(double vdWDistance, double distance) const{
+double Cndo2::GetVdWDampingValue1stDerivative(double vdWDistance, double distance) const{
    double dampingFactor = Parameters::GetInstance()->GetVdWDampingFactorSCF();
    return (dampingFactor/vdWDistance)
          *exp(-1.0*dampingFactor*(distance/vdWDistance - 1.0))
@@ -438,7 +481,7 @@ double Cndo2::GetVdwDampingValue1stDerivative(double vdWDistance, double distanc
 }
 
 // See damping function in (2) in [G_2004] ((11) in [G_2006])
-double Cndo2::GetVdwDampingValue2ndDerivative(double vdWDistance, double distance) const{
+double Cndo2::GetVdWDampingValue2ndDerivative(double vdWDistance, double distance) const{
    double dampingFactor = Parameters::GetInstance()->GetVdWDampingFactorSCF();
    double exponent = -1.0*dampingFactor*(distance/vdWDistance - 1.0);
    double pre = dampingFactor/vdWDistance;
@@ -454,7 +497,7 @@ double Cndo2::GetDiatomVdWCorrectionEnergy(const Atom& atomA, const Atom& atomB)
    double tmpSum = atomA.GetVdWCoefficient()+atomB.GetVdWCoefficient();
    if(tmpSum<=0e0){return 0e0;}
    double vdWCoefficients = 2.0*atomA.GetVdWCoefficient()*atomB.GetVdWCoefficient()/tmpSum;
-   double damping = this->GetVdwDampingValue(vdWDistance, distance);
+   double damping = this->GetVdWDampingValue(vdWDistance, distance);
    double scalingFactor = Parameters::GetInstance()->GetVdWScalingFactorSCF();
    return -1.0*scalingFactor*vdWCoefficients*damping
          /(distance*distance*distance*distance*distance*distance);
@@ -470,8 +513,8 @@ double Cndo2::GetDiatomVdWCorrection1stDerivative(const Atom& atomA, const Atom&
    if(tmpSum<=0e0){return 0e0;}
    double vdWCoefficients = 2.0*atomA.GetVdWCoefficient()*atomB.GetVdWCoefficient()/tmpSum;
    double dampingFactor = Parameters::GetInstance()->GetVdWDampingFactorSCF();
-   double damping = this->GetVdwDampingValue(vdWDistance, distance);
-   double damping1stDerivative = this->GetVdwDampingValue1stDerivative(vdWDistance, distance);
+   double damping = this->GetVdWDampingValue(vdWDistance, distance);
+   double damping1stDerivative = this->GetVdWDampingValue1stDerivative(vdWDistance, distance);
    double value=0.0;
    double tmp = distance*distance*distance*distance*distance*distance;
    value += 6.0*damping/(tmp*distance) 
@@ -498,9 +541,9 @@ double Cndo2::GetDiatomVdWCorrection2ndDerivative(const Atom& atomA,
    if(tmpSum<=0e0){return 0e0;}
    double vdWCoefficients = 2.0*atomA.GetVdWCoefficient()*atomB.GetVdWCoefficient()/tmpSum;
    double dampingFactor = Parameters::GetInstance()->GetVdWDampingFactorSCF();
-   double damping = this->GetVdwDampingValue(vdWDistance, distance);
-   double damping1stDerivative = this->GetVdwDampingValue1stDerivative(vdWDistance, distance);
-   double damping2ndDerivative = this->GetVdwDampingValue2ndDerivative(vdWDistance, distance);
+   double damping = this->GetVdWDampingValue(vdWDistance, distance);
+   double damping1stDerivative = this->GetVdWDampingValue1stDerivative(vdWDistance, distance);
+   double damping2ndDerivative = this->GetVdWDampingValue2ndDerivative(vdWDistance, distance);
 
    double dis6 = distance*distance*distance*distance*distance*distance;
    double tmp1 = -6.0*damping             /(dis6*distance) 

--- a/src/cndo/Cndo2.h
+++ b/src/cndo/Cndo2.h
@@ -63,6 +63,7 @@ protected:
    std::string errorMessageMoleculeNotSet;
    std::string errorMessageOddTotalValenceElectrions;
    std::string errorMessageNotEnebleAtomType;
+   std::string errorMessageNotEnebleAtomTypeVdW;
    std::string errorMessageCoulombInt;
    std::string errorMessageExchangeInt;
    std::string errorMessageMolecularIntegralElement;
@@ -341,6 +342,7 @@ private:
    double elecSCFEnergy;
    double bondingAdjustParameterK[2]; //see (3.79) in J. A. Pople book
    double** gammaAB;
+   std::vector<MolDS_base::AtomType> enableAtomTypesVdW;
    class ReducedOverlapAOsParameters : private MolDS_base::Uncopyable{
    public:
       // use Y[na][nb][la][lb][m][i][j] 
@@ -368,10 +370,12 @@ private:
                           double const* normalForceConstants, 
                           const MolDS_base::Molecule& molecule) const;
    void CalcCoreRepulsionEnergy();
+   void SetEnableAtomTypesVdW();
+   void CheckEnableAtomTypeVdW(const MolDS_base::Molecule& molecule) const;
    void CalcVdWCorrectionEnergy();
-   double GetVdwDampingValue(double vdWDistance, double distance) const;
-   double GetVdwDampingValue1stDerivative(double vdWDistance, double distance) const;
-   double GetVdwDampingValue2ndDerivative(double vdWDistance, double distance) const;
+   double GetVdWDampingValue(double vdWDistance, double distance) const;
+   double GetVdWDampingValue1stDerivative(double vdWDistance, double distance) const;
+   double GetVdWDampingValue2ndDerivative(double vdWDistance, double distance) const;
    void CalcElectronicDipoleMomentGroundState(double*** electronicTransitionDipoleMoments,
                                               double const* const* const* cartesianMatrix,
                                               const MolDS_base::Molecule& molecule, 

--- a/src/indo/Indo.cpp
+++ b/src/indo/Indo.cpp
@@ -72,6 +72,8 @@ void Indo::SetMessages(){
       = "Error in indo::Indo::SetMolecule: Total number of valence electrons is odd. totalNumberValenceElectrons=";
    this->errorMessageNotEnebleAtomType  
       = "Error in indo::Indo::CheckEnableAtomType: Non available atom is contained.\n";
+   this->errorMessageNotEnebleAtomTypeVdW
+      = "Error in indo::Indo::CheckEnableAtomTypeVdW: Non available atom to add VdW correction is contained.\n";
    this->errorMessageCoulombInt = "Error in base_indo::Indo::GetCoulombInt: Invalid orbitalType.\n";
    this->errorMessageExchangeInt = "Error in base_indo::Indo::GetExchangeInt: Invalid orbitalType.\n";
    this->errorMessageMolecularIntegralElement

--- a/src/mndo/Mndo.cpp
+++ b/src/mndo/Mndo.cpp
@@ -154,6 +154,8 @@ void Mndo::SetMessages(){
       = "Error in mndo::Mndo::SetMolecule: Total number of valence electrons is odd. totalNumberValenceElectrons=";
    this->errorMessageNotEnebleAtomType  
       = "Error in mndo::Mndo::CheckEnableAtomType: Non available atom is contained.\n";
+   this->errorMessageNotEnebleAtomTypeVdW
+      = "Error in mndo::Mndo::CheckEnableAtomTypeVdW: Non available atom to add VdW correction is contained.\n";
    this->errorMessageCoulombInt = "Error in base_mndo::Mndo::GetCoulombInt: Invalid orbitalType.\n";
    this->errorMessageExchangeInt = "Error in base_mndo::Mndo::GetExchangeInt: Invalid orbitalType.\n";
    this->errorMessageCalcCISMatrix

--- a/src/pm3/Pm3.cpp
+++ b/src/pm3/Pm3.cpp
@@ -79,6 +79,8 @@ void Pm3::SetMessages(){
       = "Error in pm3::Pm3::SetMolecule: Total number of valence electrons is odd. totalNumberValenceElectrons=";
    this->errorMessageNotEnebleAtomType  
       = "Error in pm3::Pm3::CheckEnableAtomType: Non available atom is contained.\n";
+   this->errorMessageNotEnebleAtomTypeVdW
+      = "Error in pm3::Pm3::CheckEnableAtomTypeVdW: Non available atom to add VdW correction is contained.\n";
    this->errorMessageCoulombInt = "Error in base_pm3::Pm3::GetCoulombInt: Invalid orbitalType.\n";
    this->errorMessageExchangeInt = "Error in base_pm3::Pm3::GetExchangeInt: Invalid orbitalType.\n";
    this->errorMessageCalcCISMatrix

--- a/src/pm3/Pm3D.cpp
+++ b/src/pm3/Pm3D.cpp
@@ -80,6 +80,8 @@ void Pm3D::SetMessages(){
       = "Error in pm3::Pm3D::SetMolecule: Total number of valence electrons is odd. totalNumberValenceElectrons=";
    this->errorMessageNotEnebleAtomType  
       = "Error in pm3::Pm3D::CheckEnableAtomType: Non available atom is contained.\n";
+   this->errorMessageNotEnebleAtomTypeVdW
+      = "Error in pm3::Pm3D::CheckEnableAtomTypeVdW: Non available atom to add VdW correction is contained.\n";
    this->errorMessageCoulombInt = "Error in pm3::Pm3D::GetCoulombInt: Invalid orbitalType.\n";
    this->errorMessageExchangeInt = "Error in pm3::Pm3D::GetExchangeInt: Invalid orbitalType.\n";
    this->errorMessageCalcCISMatrix

--- a/src/pm3/Pm3Pddg.cpp
+++ b/src/pm3/Pm3Pddg.cpp
@@ -78,6 +78,8 @@ void Pm3Pddg::SetMessages(){
       = "Error in pm3::Pm3Pddg::SetMolecule: Total number of valence electrons is odd. totalNumberValenceElectrons=";
    this->errorMessageNotEnebleAtomType  
       = "Error in pm3::Pm3Pddg::CheckEnableAtomType: Non available atom is contained.\n";
+   this->errorMessageNotEnebleAtomTypeVdW
+      = "Error in pm3::Pm3Pddg::CheckEnableAtomTypeVdW: Non available atom to add VdW correction is contained.\n";
    this->errorMessageCoulombInt = "Error in base_pm3::Pm3Pddg::GetCoulombInt: Invalid orbitalType.\n";
    this->errorMessageExchangeInt = "Error in base_pm3::Pm3Pddg::GetExchangeInt: Invalid orbitalType.\n";
    this->errorMessageCalcCISMatrix

--- a/src/zindo/ZindoS.cpp
+++ b/src/zindo/ZindoS.cpp
@@ -155,6 +155,8 @@ void ZindoS::SetMessages(){
       = "Error in zindo::ZindoS::SetMolecule: Total number of valence electrons is odd. totalNumberValenceElectrons=";
    this->errorMessageNotEnebleAtomType  
       = "Error in zindo::ZindoS::CheckEnableAtomType: Non available atom is contained.\n";
+   this->errorMessageNotEnebleAtomTypeVdW
+      = "Error in zindo::ZindoS::CheckEnableAtomTypeVdW: Non available atom to add VdW correction is contained.\n";
    this->errorMessageCoulombInt = "Error in zindo::ZindoS::GetCoulombInt: Invalid orbitalType.\n";
    this->errorMessageExchangeInt = "Error in zindo::ZindoS::GetExchangeInt: Invalid orbitalType.\n";
    this->errorMessageNishimotoMataga = "Error in zindo::ZindoS::GetNishimotoMatagaTwoEleInt: Invalid orbitalType.\n";