提交: 03a53ba3358a419e44799fb3fcec35f75b1cff6d - MolDS (git)

Commit MetaInfo

• 修订版: 03a53ba3358a419e44799fb3fcec35f75b1cff6d (tree)
• 时间: 2013-01-22 05:19:07
• 作者: Katsuhiko Nishimra <ktns.87@gmai...>
• Commiter: Katsuhiko Nishimra

Log Message

Merge trunk into automake. #28588

git-svn-id: https://svn.sourceforge.jp/svnroot/molds/branches/automake@1241 1136aad2-a195-0410-b898-f5ea1d11b9d8

更改概述

• modified: src/base/InputParser.cpp
• modified: src/base/InputParser.h
• modified: src/base/MathUtilities.cpp
• modified: src/base/MathUtilities.h
• modified: src/base/Parameters.cpp
• modified: src/base/Parameters.h
• modified: src/cndo/Cndo2.cpp
• modified: src/cndo/Cndo2.h
• modified: src/zindo/ZindoS.cpp
• modified: src/zindo/ZindoS.h

差异

--- a/src/base/InputParser.cpp
+++ b/src/base/InputParser.cpp
@@ -311,6 +311,7 @@ void InputParser::SetMessages(){
    this->stringCISAllTransitionDipoleMoments = "all_transition_dipole_moments";
    this->stringCISNumPrintCoefficients       = "num_print_coefficients";
    this->stringCISMulliken                   = "mulliken";
+   this->stringCISUnpairedPop                = "unpaired_electron_population";
 
    // Memory
    this->stringMemory          = "memory";
@@ -823,6 +824,16 @@ int InputParser::ParseConditionsCIS(vector<string>* inputTerms, int parseIndex)
          }
          parseIndex++;
       }
+      // unpaired electron population
+      if((*inputTerms)[parseIndex].compare(this->stringCISUnpairedPop) == 0){
+          if((*inputTerms)[parseIndex+1].compare(this->stringYES) == 0){
+             Parameters::GetInstance()->SetRequiresUnpairedPopCIS(true);
+          }
+          else {
+             Parameters::GetInstance()->SetRequiresUnpairedPopCIS(false);
+          }
+          parseIndex++;   
+      }
       parseIndex++;   
    }
    return parseIndex;

--- a/src/base/InputParser.h
+++ b/src/base/InputParser.h
@@ -236,6 +236,7 @@ private:
    std::string stringCISAllTransitionDipoleMoments;
    std::string stringCISNumPrintCoefficients;
    std::string stringCISMulliken;
+   std::string stringCISUnpairedPop;
    // Memory
    std::string stringMemory;
    std::string stringMemoryEnd;

--- a/src/base/MathUtilities.cpp
+++ b/src/base/MathUtilities.cpp
@@ -100,40 +100,40 @@ template <typename T> T min(T a, T b){
 }
 
 // rotating matrix
-void CalcRotatingMatrix(double matrix[][3], double sita, CartesianType cartesianType){
+void CalcRotatingMatrix(double matrix[][3], double theta, CartesianType cartesianType){
    if(cartesianType == XAxis){
       matrix[0][0] = 1.0;
       matrix[0][1] = 0.0;
       matrix[0][2] = 0.0;
 
       matrix[1][0] = 0.0;
-      matrix[1][1] = cos(sita);
-      matrix[1][2] = sin(sita);
+      matrix[1][1] = cos(theta);
+      matrix[1][2] = sin(theta);
 
       matrix[2][0] = 0.0;
-      matrix[2][1] = -sin(sita);
-      matrix[2][2] = cos(sita);
+      matrix[2][1] = -sin(theta);
+      matrix[2][2] = cos(theta);
    }
    else if(cartesianType == YAxis){
-      matrix[0][0] = cos(sita);
+      matrix[0][0] = cos(theta);
       matrix[0][1] = 0.0;
-      matrix[0][2] = -sin(sita);
+      matrix[0][2] = -sin(theta);
 
       matrix[1][0] = 0.0;
       matrix[1][1] = 1.0;
       matrix[1][2] = 0.0;
 
-      matrix[2][0] = sin(sita);
+      matrix[2][0] = sin(theta);
       matrix[2][1] = 0.0;
-      matrix[2][2] = cos(sita);
+      matrix[2][2] = cos(theta);
    }
    else if(cartesianType == ZAxis){
-      matrix[0][0] = cos(sita);
-      matrix[0][1] = sin(sita);
+      matrix[0][0] = cos(theta);
+      matrix[0][1] = sin(theta);
       matrix[0][2] = 0.0;
 
-      matrix[1][0] = -sin(sita);
-      matrix[1][1] = cos(sita);
+      matrix[1][0] = -sin(theta);
+      matrix[1][1] = cos(theta);
       matrix[1][2] = 0.0;
 
       matrix[2][0] = 0.0;

--- a/src/base/MathUtilities.h
+++ b/src/base/MathUtilities.h
@@ -28,7 +28,7 @@ template <typename T> T Max(T a, T b);
 // min
 template <typename T> T min(T a, T b);
 // rotating matrix
-void CalcRotatingMatrix(double matrix[][3], double sita, CartesianType cartesianType);
+void CalcRotatingMatrix(double matrix[][3], double theta, CartesianType cartesianType);
 // calculate determinant of the matrix. Note taht the matrix will be destroid
 double GetDeterminant(double** matrix, int dim);
 }

--- a/src/base/Parameters.cpp
+++ b/src/base/Parameters.cpp
@@ -1,6 +1,7 @@
 //************************************************************************//
 // Copyright (C) 2011-2012 Mikiya Fujii                                   // 
 // Copyright (C) 2012-2012 Katsuhiko Nishimra                             // 
+// Copyright (C) 2012-2013 Michihiro Okuyama
 //                                                                        // 
 // This file is part of MolDS.                                            // 
 //                                                                        // 
@@ -171,6 +172,7 @@ void Parameters::SetDefaultValues(){
    this->numberPrintCoefficientsCIS            = 1;
    this->requiresExcitonEnergiesCIS            = false;
    this->requiresAllTransitionDipoleMomentsCIS = false;
+   this->requiresUnpairedPopCIS                = false;
    // Memory
    this->limitHeapMemory = 256;
    // MD
@@ -705,6 +707,14 @@ bool Parameters::RequiresMullikenCIS() const{
            0<this->electronicStateIndecesMullikenCIS->size());
 }
 
+void Parameters::SetRequiresUnpairedPopCIS(bool requiresUnpairedPopCIS){
+   this->requiresUnpairedPopCIS = requiresUnpairedPopCIS;
+}
+
+bool Parameters::RequiresUnpairedPopCIS() const{
+   return this->requiresUnpairedPopCIS; 
+}
+
 // methods for Memory
 double Parameters::GetLimitHeapMemory() const{
    return this->limitHeapMemory;

--- a/src/base/Parameters.h
+++ b/src/base/Parameters.h
@@ -1,6 +1,7 @@
 //************************************************************************//
 // Copyright (C) 2011-2012 Mikiya Fujii                                   // 
 // Copyright (C) 2012-2012 Katsuhiko Nishimra                             // 
+// Copyright (C) 2012-2013 Michihiro Okuyama
 //                                                                        // 
 // This file is part of MolDS.                                            // 
 //                                                                        // 
@@ -141,6 +142,8 @@ public:
    std::vector<int>* GetElectronicStateIndecesMullikenCIS() const;
    void              AddElectronicStateIndexMullikenCIS(int electronicStateIndex);
    bool              RequiresMullikenCIS() const;
+   void              SetRequiresUnpairedPopCIS(bool requiresUnpairedPopCIS); 
+   bool              RequiresUnpairedPopCIS() const; 
    // Memory
    double GetLimitHeapMemory() const;
    void   SetLimitHeapMemory(double limitHeap);
@@ -285,6 +288,7 @@ private:
    bool              requiresExcitonEnergiesCIS;
    bool              requiresAllTransitionDipoleMomentsCIS;
    std::vector<int>* electronicStateIndecesMullikenCIS;
+   bool              requiresUnpairedPopCIS;
    // Memory
    double limitHeapMemory;
    // MD

--- a/src/cndo/Cndo2.cpp
+++ b/src/cndo/Cndo2.cpp
@@ -1,5 +1,6 @@
 //************************************************************************//
 // Copyright (C) 2011-2012 Mikiya Fujii                                   // 
+// Copyright (C) 2012-2013 Michihiro Okuyama
 //                                                                        // 
 // This file is part of MolDS.                                            // 
 //                                                                        // 
@@ -73,6 +74,7 @@ Cndo2::Cndo2(){
    this->orbitalElectronPopulationCIS = NULL;
    this->atomicElectronPopulation     = NULL;
    this->atomicElectronPopulationCIS  = NULL;
+   this->atomicUnpairedPopulationCIS  = NULL;
    this->overlapAOs = NULL;
    this->twoElecTwoCore = NULL;
    this->cartesianMatrix = NULL;
@@ -206,6 +208,8 @@ void Cndo2::SetMessages(){
    this->messageMullikenAtomsSCF   = "\tMulliken charge(SCF):";
    this->messageMullikenAtoms      = "\tMulliken charge:";
    this->messageMullikenAtomsTitle = "\t\t\t\t| k-th eigenstate | i-th atom | atom type | core charge[a.u.] | Mulliken charge[a.u.]| \n";
+   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";

--- a/src/cndo/Cndo2.h
+++ b/src/cndo/Cndo2.h
@@ -1,5 +1,6 @@
 //************************************************************************//
 // Copyright (C) 2011-2012 Mikiya Fujii                                   // 
+// Copyright (C) 2012-2013 Michihiro Okuyama
 //                                                                        // 
 // This file is part of MolDS.                                            // 
 //                                                                        // 
@@ -86,6 +87,8 @@ protected:
    std::string messageOmpElapsedTimeSCF;
    std::string messageMullikenAtoms;
    std::string messageMullikenAtomsTitle;
+   std::string messageUnpairedAtoms;
+   std::string messageUnpairedAtomsTitle;
    std::string messageUnitSec; 
    std::vector<MolDS_base::AtomType> enableAtomTypes;
    MolDS_base::Molecule* molecule;
@@ -99,6 +102,7 @@ protected:
    double***    orbitalElectronPopulationCIS; 
    double*      atomicElectronPopulation; //P_{AB} of (3.21) in J. A. Pople book.
    double**     atomicElectronPopulationCIS; 
+   double**     atomicUnpairedPopulationCIS; 
    double**     overlapAOs; // overlap integral between AOs
    double****** twoElecTwoCore;
    double***    cartesianMatrix; // cartesian matrix represented by AOs

--- a/src/zindo/ZindoS.cpp
+++ b/src/zindo/ZindoS.cpp
@@ -1,5 +1,6 @@
 //************************************************************************//
 // Copyright (C) 2011-2012 Mikiya Fujii                                   // 
+// Copyright (C) 2012-2013 Michihiro Okuyama
 //                                                                        // 
 // This file is part of MolDS.                                            // 
 //                                                                        // 
@@ -1021,6 +1022,10 @@ void ZindoS::CalcCISProperties(){
    this->CalcAtomicElectronPopulationCIS(&this->atomicElectronPopulationCIS,
                                          this->orbitalElectronPopulationCIS, 
                                          *this->molecule);
+   // atomic unpaired electron population
+   this->CalcAtomicUnpairedPopulationCIS(&this->atomicUnpairedPopulationCIS,
+                                         this->orbitalElectronPopulationCIS, 
+                                         *this->molecule);
 }
 
 void ZindoS::CalcElectronicDipoleMomentsExcitedState(double*** electronicTransitionDipoleMoments,
@@ -1407,6 +1412,58 @@ void ZindoS::CalcAtomicElectronPopulationCIS(double*** atomicElectronPopulationC
    }
 }
 
+void ZindoS::CalcAtomicUnpairedPopulationCIS(double*** atomicUnpairedPopulationCIS,
+                                             double const* const* const* orbitalElectronPopulationCIS, 
+                                             const Molecule& molecule) const{
+   if(!Parameters::GetInstance()->RequiresMullikenCIS()){
+      return;
+   }
+   if(!Parameters::GetInstance()->RequiresUnpairedPopCIS()){
+      return;
+   }
+   int totalNumberAtoms = molecule.GetNumberAtoms();
+   vector<int>* elecStates = Parameters::GetInstance()->GetElectronicStateIndecesMullikenCIS();
+   // malloc or initialize free exciton energies
+   if(*atomicUnpairedPopulationCIS == NULL){
+      MallocerFreer::GetInstance()->Malloc<double>(atomicUnpairedPopulationCIS, 
+                                                   elecStates->size(),
+                                                   totalNumberAtoms);
+   }
+   else{
+      MallocerFreer::GetInstance()->Initialize<double>(*atomicUnpairedPopulationCIS,
+                                                       elecStates->size(),
+                                                       totalNumberAtoms);
+   }
+   // calc atomic electron population
+   for(int k=0; k<elecStates->size(); k++){
+      stringstream ompErrors;
+#pragma omp parallel for schedule(auto)
+      for(int a=0; a<totalNumberAtoms; a++){
+         try{
+            int firstAOIndex = molecule.GetAtom(a)->GetFirstAOIndex();
+            int numberAOs = molecule.GetAtom(a)->GetValenceSize();
+            (*atomicUnpairedPopulationCIS)[k][a] = 0.0; 
+            for(int i=firstAOIndex; i<firstAOIndex+numberAOs; i++){
+               double orbitalSquarePopulation = 0.0; 
+               int    totalNumberAOs = molecule.GetTotalNumberAOs(); 
+               for(int j=0; j<totalNumberAOs; j++) {
+                  orbitalSquarePopulation += orbitalElectronPopulationCIS[k][i][j] * orbitalElectronPopulationCIS[k][j][i]; 
+               }
+               (*atomicUnpairedPopulationCIS)[k][a] += 2.0 * orbitalElectronPopulationCIS[k][i][i] - orbitalSquarePopulation;
+            }
+         }
+         catch(MolDSException ex){
+#pragma omp critical
+            ompErrors << ex.what() << endl ;
+         }
+      }
+      // Exception throwing for omp-region
+      if(!ompErrors.str().empty()){
+         throw MolDSException(ompErrors.str());
+      }
+   }
+}
+
 void ZindoS::OutputCISResults() const{
    int numberActiveOcc = Parameters::GetInstance()->GetActiveOccCIS();
    int numberActiveVir = Parameters::GetInstance()->GetActiveVirCIS();
@@ -1442,6 +1499,9 @@ void ZindoS::OutputCISResults() const{
    // output mulliken population
    this->OutputCISMulliken();
 
+   // output unpaired electron population
+   this->OutputCISUnpairedPop(); 
+
    // output exciton energies
    if(Parameters::GetInstance()->RequiresExcitonEnergiesCIS()){
       this->OutputLog(this->messageExcitonEnergiesCIS);
@@ -1589,6 +1649,29 @@ void ZindoS::OutputCISMulliken() const{
    }
 }
 
+void ZindoS::OutputCISUnpairedPop() const{
+   if(!Parameters::GetInstance()->RequiresMullikenCIS()){
+      return;
+   }
+   if(!Parameters::GetInstance()->RequiresUnpairedPopCIS()){
+      return;
+   }
+   int totalNumberAtoms = this->molecule->GetNumberAtoms();
+   this->OutputLog(this->messageUnpairedAtomsTitle);
+   vector<int>* elecStates = Parameters::GetInstance()->GetElectronicStateIndecesMullikenCIS();
+   for(int k=0; k<elecStates->size(); k++){
+      for(int a=0; a<totalNumberAtoms; a++){
+         const Atom& atom = *this->molecule->GetAtom(a);
+         this->OutputLog(boost::format("%s\t%d\t%d\t%s\t%e\n") % this->messageUnpairedAtoms
+                                                               % (*elecStates)[k]
+                                                               % a
+                                                               % AtomTypeStr(atom.GetAtomType())
+                                                               % this->atomicUnpairedPopulationCIS[k][a]);
+      }
+      this->OutputLog("\n");
+   }
+}
+
 void ZindoS::SortCISEigenVectorCoefficients(vector<CISEigenVectorCoefficient>* cisEigenVectorCoefficients,
                                             double* cisEigenVector) const{
    for(int l=0; l<this->matrixCISdimension; l++){

--- a/src/zindo/ZindoS.h
+++ b/src/zindo/ZindoS.h
@@ -1,5 +1,6 @@
 //************************************************************************//
 // Copyright (C) 2011-2012 Mikiya Fujii                                   // 
+// Copyright (C) 2012-2013 Michihiro Okuyama
 //                                                                        // 
 // This file is part of MolDS.                                            // 
 //                                                                        // 
@@ -139,6 +140,7 @@ private:
    void OutputCISDipole() const;
    void OutputCISTransitionDipole() const;
    void OutputCISMulliken() const;
+   void OutputCISUnpairedPop() const;
    void CalcFreeExcitonEnergies(double** freeExcitonEnergiesCIS, 
                                 const MolDS_base::Molecule& molecule, 
                                 double const* energiesMO, 
@@ -152,6 +154,9 @@ private:
    void CalcAtomicElectronPopulationCIS(double*** atomicElectronPopulationCIS,
                                         double const* const* const* orbitalElectronPopulationCIS, 
                                         const MolDS_base::Molecule& molecule) const;
+   void CalcAtomicUnpairedPopulationCIS(double*** atomicUnpairedPopulationCIS,
+                                        double const* const* const* orbitalElectronPopulationCIS, 
+                                        const MolDS_base::Molecule& molecule) const; 
    void CalcElectronicDipoleMomentsExcitedState(double*** electronicTransitionDipoleMoments,
                                                 double const* const* fockMatrix,
                                                 double const* const* matrixCIS,