00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #include "QCCommon.hpp"
00020 #include "QCReader.hpp"
00021 #include "QCMndo.hpp"
00022 #include "QCAm1.hpp"
00023 #include "QCPm3.hpp"
00024 #include "QCMndoParam.hpp"
00025 #include "QCAm1Param.hpp"
00026 #include "QCPm3Param.hpp"
00027 #include "QCSCFAlgo.hpp"
00028 #include "QCManager.hpp"
00029 #include "QCGlobalSystem.hpp"
00030 #include "QCMemory.hpp"
00031 #include "QCSymMatrix.hpp"
00032
00033 #include "QCTools.hpp"
00034
00035
00036
00037
00038 #ifndef QC_NO_DEBUG
00039 # define VERBOSE_DERIV
00040 # define DEBUG_DERIV
00041 #endif
00042
00043 template <class TPSolver>
00044 bool QCGlobalSystem<TPSolver>::checkData(const QCGeneralData& data){
00045 QC_TRACE_INIT("BEGIN QCGlobalSystem<TPSolver>::checkdata(const QCGeneralData& data)");
00046
00047 QC_TRACE_INIT("END QCGlobalSystem<TPSolver>::checkdata(const QCGeneralData& data)");
00048 return true ;
00049 }
00053 template <class TPSolver>
00054 template <class TPManager>
00055 void
00056 QCGlobalSystem<TPSolver>::computeEnergyDerivatives (TPManager& manager) {
00057
00058 #ifdef VERBOSE_DERIV
00059 cout << endl << " ** COMPUTE ENERGY GRADIENT ** " << endl;
00060 #endif
00061
00062 QCMemory& QCRestrict memory = manager.getMemory();
00063
00064 const QCGeneralData& data = manager.getGeneralData();
00065
00066 typedef typename TPManager::TModel TModel;
00067 TModel& QCRestrict model = manager.getModel();
00068
00069 QCFloat deltaQ = data.getDeltaQ();
00070
00071 QCFloat * QCRestrict dRepIntegralSSVector = NULL;
00072
00073
00074 typename TPManager::TSystem::QCIterator systemIter = begin();
00075
00076 QCPoint3D shift;
00077
00078 dRepIntegralSSVector = memory.template takeTmpMem<QCFloat>(nbAtoms*(nbAtoms-1)/2);
00079 memset(dRepIntegralSSVector, 0, (nbAtoms*(nbAtoms-1))/2*sizeof(QCFloat));
00080 overlapS->getMatrix().takeTmpMem(memory);
00081 dOverlapSOnDr->getMatrix().takeTmpMem(memory);
00082
00083
00084
00085
00086
00087 #ifdef DEBUG_DERIV
00088 displayDisps(cout);
00089 #endif
00090
00091
00092 shift[COORDX] = deltaQ;
00093 shift[COORDY] = QC_ZERO;
00094 shift[COORDZ] = QC_ZERO;
00095 dOverlapSOnDr->computeDerivatives(manager, systemIter, shift);
00096 model.deriveElecEnergy(manager,systemIter, shift, dRepIntegralSSVector);
00097
00098
00099
00100 shift[COORDX] = QC_ZERO;
00101 shift[COORDY] = deltaQ;
00102 shift[COORDZ] = QC_ZERO;
00103 dOverlapSOnDr->computeDerivatives(manager, systemIter, shift);
00104 model.deriveElecEnergy(manager,systemIter,shift,dRepIntegralSSVector);
00105
00106
00107
00108 shift[COORDX] = QC_ZERO;
00109 shift[COORDY] = QC_ZERO;
00110 shift[COORDZ] = deltaQ;
00111 dOverlapSOnDr->computeDerivatives(manager, systemIter, shift);
00112 model.deriveElecEnergy(manager, systemIter,shift, dRepIntegralSSVector);
00113
00114 dOverlapSOnDr->getMatrix().giveBackTmpMem(memory);
00115 overlapS->getMatrix().giveBackTmpMem(memory);
00116
00117
00118 divDisps(QC_TWO*deltaQ);
00119
00120 model.deriveCoreEnergy(manager,systemIter,dRepIntegralSSVector);
00121
00122 memory.template giveBackTmpMem<QCFloat>(dRepIntegralSSVector);
00123
00124 #ifdef DEBUG_DERIV
00125 displayDisps(cout);
00126 #endif
00127 }
00128
00129 template <class TPSolver>
00130 void
00131 QCGlobalSystem<TPSolver>::writeDensityOnFileAscii(const QCFiles & files, const string& path){
00132 QC_TRACE_OUT("BEGIN QCGlobalSystem<TPSolver>::writeDensityOnFileAscii ");
00133
00134 std::ofstream out ;
00135 std::string fileName;
00136
00137 fileName = path + "/" + files.getResultFile() + "-density.ascii";
00138
00139 out.open(fileName.c_str());
00140
00141 out.setf(std::ios_base::scientific) ; out.precision(20) ;
00142
00143 QCSymMatrix& QCRestrict D = QCModelMatrices::densityP.getMatrix();
00144 out << "BASIC "<< D.getDim() << " "<< D.getNbElems() << std::endl;
00145 for (int i = 0; i < D.getDim(); ++i) {
00146 for (int j = 0; j <=i; ++j) {
00147 out << " "<< i << " "<< j << " "<< D[i][j] << std::endl;
00148 }
00149 }
00150 out.close();
00151 QC_TRACE_OUT("END QCGlobalSystem<TPSolver>::writeDensityOnFileAscii " ) ;
00152 }
00153 template <class TPSolver>
00154 void
00155 QCGlobalSystem<TPSolver>::writeDensityOnFileBinary(const QCFiles & files, const string& path){
00156 QC_TRACE_OUT("BEGIN QCGlobalSystem<TPSolver>::writeDensityOnFileBinary ");
00157 std::ofstream out ;
00158 std::string fileName,type("SEQ ");
00159 fileName = path + "/" + files.getResultFile() + "-density.bin";
00160 out.open(fileName.c_str(), ios::out | ios::binary);
00161 int val[2];
00162
00163
00164 QCSymMatrix& QCRestrict D = QCModelMatrices::densityP.getMatrix();
00165 val[0] = D.getDim(); val[1] = D.getNbElems()*sizeof(QCFloat) ;
00166 out.write (type.c_str(), 4);
00167 out.write ((char*)(&val[0]), 2*sizeof(int));
00168 out.write ((char*)(D.getElems()), val[1]);
00169
00170 out.close();
00171 QC_TRACE_OUT("END QCGlobalSystem<TPSolver>::writeDensityOnFileBinary " ) ;
00172 }
00173 template <class TPSolver>
00174 void
00175 QCGlobalSystem<TPSolver>::readDensityFromFile(const QCFiles & files, const string& path){
00176
00177
00178 QC_TRACE_INIT("BEGIN QCGlobalSystem<TPSolver>::readDensityFromFile ");
00179 std::string fileName;
00180 fileName = files.getDensityFile();
00181
00182
00183 if(fileName.find(".bin")!= std::string::npos){
00184 this->readDensityFromFileBinary( files, path);
00185 }
00186 else if(fileName.find(".ascii",0)!= std::string::npos){
00187 readGlobalDensityFromAsciiFile(fileName,QCModelMatrices::densityP.getMatrix());
00188 }
00189 else {
00190 std::cerr <<" Wrong file for the density matrix. the extension must be .bin or .ascii"<<std::endl;
00191 exit(EXIT_FAILURE) ;
00192 }
00193 QC_TRACE_INIT("END QCGlobalSystem<TPSolver>::readDensityFromFile ");
00194 }
00195
00196 template <class TPSolver>
00197 void
00198 QCGlobalSystem<TPSolver>::readDensityFromFileAscii(const QCFiles & files, const string& path){
00199 QC_TRACE_INIT("BEGIN QCGlobalSystem<TPSolver>::readDensityFromFileAscii ") ;
00200 std::ifstream data ;
00201 std::string fileName, type;
00202 int dim = 0, nbelement = 0;
00203 fileName = files.getDensityFile();
00204 data.open(fileName.c_str());
00205 if (!data){
00206 std::cerr << "Error to open density Matrix in File " << fileName <<std::endl;
00207 exit(EXIT_FAILURE) ;
00208 }
00209
00210 QCSymMatrix& QCRestrict D = QCModelMatrices::densityP.getMatrix();
00211
00212 data >> type >> dim >> nbelement;
00213 if(type != "SEQ") {
00214 std::cerr << "Bad type of file " << type << " must be seq FOR binary file "<<std::endl;
00215 exit (EXIT_FAILURE) ;
00216 }
00217 if ( dim != D.getDim() ){
00218 std::cerr << "Bad dimension of the density Matrix in File " << fileName <<std::endl
00219 << " Dimension in File : " << dim <<" Dimension given by the model" << D.getDim() <<std::endl ;
00220 exit (EXIT_FAILURE) ;
00221 }
00222 int i, j ;
00223 QCFloat val ;
00224 for (int k = 0; k < nbelement ; ++k) {
00225 data >> i >> j >>val ;
00226 if (i < j ){
00227 std::cerr <<" Error we must read a Symetric Matrix " <<std::endl ;
00228 exit(EXIT_FAILURE) ;
00229 }
00230 D[i][j] = val ;
00231 }
00232 data.close();
00233 QC_TRACE_INIT("END QCGlobalSystem<TPSolver>::readDensityFromFileAscii ")
00234 }
00235 template <class TPSolver>
00236 void
00237 QCGlobalSystem<TPSolver>::readDensityFromFileBinary(const QCFiles & files, const string& path){
00238 QC_TRACE_INIT("BEGIN QCGlobalSystem<TPSolver>::readDensityFromFileBinary ") ;
00239 std::ifstream data ;
00240 std::string fileName, tt ;
00241 char type[4];
00242 fileName = files.getDensityFile();
00243 data.open(fileName.c_str(), ios::in|ios::binary);
00244 if (!data){
00245 std::cerr << "Error to open density Matrix in File " << fileName <<std::endl;
00246 exit(EXIT_FAILURE) ;
00247 }
00248
00249 QCSymMatrix& QCRestrict D = QCModelMatrices::densityP.getMatrix();
00250 int val[2] ;
00251
00252 data.read(type, 4); tt = type ;
00253 data.read((char*)(&val[0]), 2*sizeof(int));
00254 if(tt != "SEQ") {
00255 std::cerr << "Bad type of file " << type << " must be seq FOR binary file "<<std::endl;
00256 exit (EXIT_FAILURE) ;
00257 }
00258 if ( val[0] != D.getDim() ){
00259 std::cerr << "Bad dimension of the density Matrix in File " << fileName <<std::endl
00260 << " Dimension in File : " << val[0] <<" Dimension given by the model" << D.getDim() <<std::endl ;
00261 exit (EXIT_FAILURE) ;
00262 }
00263 data.read((char*)(D.getElems()), val[1]);
00264
00265 data.close();
00266 QC_TRACE_INIT("END QCGlobalSystem<TPSolver>::readDensityFromFileBinary ")
00267 }
00268 template <class TPSolver>
00269 void
00270 QCGlobalSystem<TPSolver>::writeSystemOnFile(const QCFiles & files, const string& path){
00271 QC_TRACE_OUT("BEGIN QCGlobalSystem<TPSolver>::writeSystemOnFile ") ;
00272 std::ofstream out ;
00273 std::string fileName;
00274
00275 fileName = path + "/" + files.getResultFile() + ".system";
00276 out.open(fileName.c_str());
00277
00278 QCSystem::writeQCIAtomicSystem(out) ;
00279 QC_TRACE_OUT("END QCGlobalSystem<TPSolver>::writeSystemOnFile ") ;
00280 }
00284 template class QCGlobalSystem<QCSCFAlgo>;
00285
00286
00287 QCMANAGER_METH_EXPL_INST_GLOBAL(void QCGlobalSystem<QCSCFAlgo>::computeEnergyDerivatives);