Add a sample code for coupled MC.

git-svn-id: https://subversion.assembla.com/svn/qmcdev/trunk@5918 e5b18d87-469d-4833-9cc0-8cdfa06e9491

src/QMCApp/CMakeLists.txt

@@ -11,6 +11,7 @@ SET(QMCAPPDIR
     QMCAppBase.cpp
     QMCDriverFactory.cpp
     QMCMain.cpp
+    CoupledMC.cpp
 )
 ADD_LIBRARY(qmc ${QMCAPPDIR})
 

src/QMCApp/CoupledMC.cpp

@@ -0,0 +1,109 @@
+//////////////////////////////////////////////////////////////////
+// (c) Copyright 2003- by Jeongnim Kim
+//////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////
+// -*- C++ -*-
+/**@file QMCMain.cpp
+ * @brief Implments QMCMain operators.
+ */
+#include "QMCApp/QMCMain.h"
+#include "QMCApp/ParticleSetPool.h"
+#include "QMCApp/WaveFunctionPool.h"
+#include "QMCApp/HamiltonianPool.h"
+#include "QMCWaveFunctions/TrialWaveFunction.h"
+#include "QMCHamiltonians/QMCHamiltonian.h"
+#include "ParticleBase/RandomSeqGenerator.h"
+#include "Utilities/OhmmsInfo.h"
+#include "Utilities/Timer.h"
+#include "QMCDrivers/QMCDriver.h"
+#include "QMCDrivers/CMC/DummyIonMove.h"
+#include "Message/Communicate.h"
+#include "Message/OpenMP.h"
+#include "qmc_common.h"
+using namespace std;
+#include "OhmmsData/AttributeSet.h"
+
+namespace qmcplusplus
+{
+
+bool QMCMain::executeCMCSection(xmlNodePtr cur)
+{
+  bool success=true;
+  string target("ion0");
+  OhmmsAttributeSet a;
+  a.add(target,"target");
+  a.put(cur);
+
+  MCWalkerConfiguration *ions = ptclPool->getWalkerSet(target);
+  TrialWaveFunction* primaryPsi=psiPool->getPrimary();
+  QMCHamiltonian* primaryH=hamPool->getPrimary();
+
+  app_log() << "QMCMain::executeCMCSection moving " << target << " by dummy move." << endl;
+  //DummyIonMove dummy(*ions,*primaryPsi,*primaryH,*hamPool,*psiPool,qmcDriver);
+  //dummy.run();
+
+  int nat = ions->getTotalNum();
+  ParticleSet::ParticlePos_t deltaR(nat);
+
+  makeGaussRandomWithEngine(deltaR,Random); //generate random displaement
+
+  //update QMC system
+  qmcSystem->update();
+
+  double logpsi1 = primaryPsi->evaluateLog(*qmcSystem);
+  cout << "logpsi1 " << logpsi1 << endl;
+
+  double eloc1  = primaryH->evaluate(*qmcSystem);
+  cout << "Local Energy " << eloc1 << endl;
+
+  for (int i=0; i<primaryH->sizeOfObservables(); i++)
+    app_log() << "  HamTest " << primaryH->getObservableName(i) << " " << primaryH->getObservable(i) << endl;
+
+  for(int iat=0; iat<nat; ++iat)
+  {
+    ions->R[iat]+=deltaR[iat];
+
+    ions->update(); //update position and distance table of itself 
+    primaryH->update_source(*ions);
+
+    qmcSystem->update();
+    double logpsi2 = primaryPsi->evaluateLog(*qmcSystem);
+    double eloc2  = primaryH->evaluate(*qmcSystem);
+
+    cout << "\nION " << iat << " " << ions->R[iat] << endl;
+    cout << "logpsi " << logpsi2 << endl;
+    cout << "Local Energy " << eloc2 << endl;
+    for (int i=0; i<primaryH->sizeOfObservables(); i++)
+      app_log() << "  HamTest " << primaryH->getObservableName(i) << " " << primaryH->getObservable(i) << endl;
+
+    ions->R[iat]-=deltaR[iat];
+    ions->update(); //update position and distance table of itself 
+    primaryH->update_source(*ions);
+
+    qmcSystem->update();
+    double logpsi3 = primaryPsi->evaluateLog(*qmcSystem);
+    double eloc3  = primaryH->evaluate(*qmcSystem);
+
+    if(abs(eloc1-eloc3)>1e-12)
+    {
+      cout << "ERROR Energies are different " << endl;
+    }
+  }
+
+  //
+  //for(int i=0; i<ions->getTotalNum(); ++i)
+  //{
+  //  ions->R[i]+=deltaR(i);
+  //  ions->update();
+  //  primaryH->update_source(*ions);
+  //  qmcDriver->run();
+  //}
+  return success;
+}
+
+}
+/***********************************************************************
+ * $RCSfilMCMain.cpp,v $   $Author: jnkim $
+ * $Revision: 5845 $   $Date: 2013-05-14 16:10:18 -0400 (Tue, 14 May 2013) $
+ * $Id: QMCMain.cpp 5845 2013-05-14 20:10:18Z jnkim $
+ ***************************************************************************/

src/QMCApp/ParticleSetPool.cpp

@@ -41,6 +41,16 @@ ParticleSetPool::ParticleSetPool(Communicate* c, const char* aname)
   myName=aname;
 }
 
+void ParticleSetPool::make_clones(int n)
+{
+  PoolType::const_iterator it(myPool.begin()), it_end(myPool.end());
+  while(it != it_end)
+  {
+    (*it).second->make_clones(n);
+    ++it;
+  }
+}
+
 ParticleSet* ParticleSetPool::getParticleSet(const string& pname)
 {
   map<string,ParticleSet*>::iterator pit(myPool.find(pname));

src/QMCApp/ParticleSetPool.h

@@ -98,6 +98,10 @@ public:
    */
   void randomize();
 
+  /** make clones for the ParticleSets of this pool
+   *    */
+  void make_clones(int n);
+
 private:
   /** global SimulationCell
    *

src/QMCApp/QMCMain.cpp

@@ -115,6 +115,10 @@ bool QMCMain::execute()
       executeLoop(cur);
       qmc_common.qmc_counter=0;
     }
+    else if(cname == "cmc")
+    {
+      executeCMCSection(cur);
+    }
   }
   m_qmcaction.clear();
   app_log() << "  Total Execution time = " << t1.elapsed() << " secs" << endl;
@@ -280,50 +284,44 @@ bool QMCMain::validateXML()
     {
       putCommunicator(cur);
     }
-    else
-      if(cname == "particleset")
+    else if(cname == "particleset")
+    {
+      ptclPool->put(cur);
+    }
+    else if(cname == "wavefunction")
+    {
+      psiPool->put(cur);
+    }
+    else if(cname == "hamiltonian")
+    {
+      hamPool->put(cur);
+    }
+    else if(cname == "include")
+    {
+      //file is provided
+      const xmlChar* a=xmlGetProp(cur,(const xmlChar*)"href");
+      if(a)
       {
-        ptclPool->put(cur);
+        pushDocument((const char*)a);
+        inputnode = processPWH(XmlDocStack.top()->getRoot());
+        popDocument();
       }
-      else
-        if(cname == "wavefunction")
-        {
-          psiPool->put(cur);
-        }
-        else
-          if(cname == "hamiltonian")
-          {
-            hamPool->put(cur);
-          }
-          else
-            if(cname == "include")
-            {
-              //file is provided
-              const xmlChar* a=xmlGetProp(cur,(const xmlChar*)"href");
-              if(a)
-              {
-                pushDocument((const char*)a);
-                inputnode = processPWH(XmlDocStack.top()->getRoot());
-                popDocument();
-              }
-            }
-            else
-              if(cname == "qmcsystem")
-              {
-                processPWH(cur);
-              }
-              else
-                if(cname == "init")
-                {
-                  InitMolecularSystem moinit(ptclPool);
-                  moinit.put(cur);
-                }
-                else
-                {
-                  //everything else goes to m_qmcaction
-                  m_qmcaction.push_back(pair<xmlNodePtr,bool>(cur,true));
-                  inputnode=false;
-                }
+    }
+    else if(cname == "qmcsystem")
+    {
+      processPWH(cur);
+    }
+    else if(cname == "init")
+    {
+      InitMolecularSystem moinit(ptclPool);
+      moinit.put(cur);
+    }
+    else
+    {
+      //everything else goes to m_qmcaction
+      m_qmcaction.push_back(pair<xmlNodePtr,bool>(cur,true));
+      inputnode=false;
+    }
     if(inputnode)
       lastInputNode=cur;
     cur=cur->next;
@@ -373,28 +371,25 @@ bool QMCMain::processPWH(xmlNodePtr cur)
     {
       ptclPool->putLattice(cur);
     }
+    else if(cname == "particleset")
+    {
+      ptclPool->putTileMatrix(cur_root);
+      ptclPool->put(cur);
+    }
+    else if(cname == "wavefunction")
+    {
+      psiPool->put(cur);
+    }
+    else if(cname == "hamiltonian")
+    {
+      hamPool->put(cur);
+    }
     else
-      if(cname == "particleset")
-      {
-        ptclPool->putTileMatrix(cur_root);
-        ptclPool->put(cur);
-      }
-      else
-        if(cname == "wavefunction")
-        {
-          psiPool->put(cur);
-        }
-        else
-          if(cname == "hamiltonian")
-          {
-            hamPool->put(cur);
-          }
-          else
-            //add to m_qmcaction
-          {
-            inputnode=false;
-            m_qmcaction.push_back(pair<xmlNodePtr,bool>(xmlCopyNode(cur,1),false));
-          }
+      //add to m_qmcaction
+    {
+      inputnode=false;
+      m_qmcaction.push_back(pair<xmlNodePtr,bool>(xmlCopyNode(cur,1),false));
+    }
     cur=cur->next;
   }
   //flush

src/QMCApp/QMCMain.h

@@ -86,6 +86,9 @@ private:
    * @return true, if a section is successfully executed.
    */
   bool executeQMCSection(xmlNodePtr cur, bool noloop=true);
+  ///execute <cmc/> element
+  bool executeCMCSection(xmlNodePtr cur);
+
 };
 }
 #endif