1patch -u -l -b -F 5 -N --suffix=.preplumed "./Makefile" << \EOF_EOF
2--- ./Makefile.preplumed
3+++ ./Makefile
4@@ -438,10 +438,15 @@
5 CXXFLAGS = $(CXXBASEFLAGS) $(CXXOPTS)
6 CXXMICFLAGS = $(CXXBASEFLAGS) $(CXXOPTS) $(CXXMICOPTS)
7 CXXTHREADFLAGS = $(CXXBASEFLAGS) $(CXXTHREADOPTS)
8 CXXSIMPARAMFLAGS = $(CXXBASEFLAGS) $(CXXSIMPARAMOPTS)
9 CXXNOALIASFLAGS = $(CXXBASEFLAGS) $(CXXNOALIASOPTS)
10+include .rootdir/Plumed.inc
11+PLUMED_PREFIX=-Wl,
12+PLUMED_LOAD1=$(patsubst "%.so",$(PLUMED_PREFIX)"%.so",$(PLUMED_LOAD))
13+PLUMED_LOAD2=$(patsubst "%.dylib",$(PLUMED_PREFIX)"%.dylib",$(PLUMED_LOAD1))
14+
15 COLVARSCXXFLAGS = $(CXXBASEFLAGS) $(CXXCOLVAROPTS) $(COPTI)$(LEPTONINCDIR) -DLEPTON -DLEPTON_USE_STATIC_LIBRARIES
16 GXXFLAGS = $(CXXBASEFLAGS) -DNO_STRSTREAM_H
17 CFLAGS = $(COPTI)$(SRCDIR) $(TCL) $(COPTS) $(RELEASE) $(EXTRADEFINES) $(TRACEOBJDEF)
18 PLUGINGCCFLAGS = $(COPTI)$(PLUGINSRCDIR) $(COPTI)$(PLUGININCDIR) $(COPTD)STATIC_PLUGIN
19 PLUGINCFLAGS = $(PLUGINGCCFLAGS) $(COPTS)
20@@ -481,10 +486,11 @@
21 	$(CUDALIB) \
22 	$(DPMTALIB) \
23 	$(DPMELIB) \
24 	$(FMMLIB) \
25 	$(TCLLIB) \
26+        $(PLUMED_LOAD2) \
27 	$(PYTHONLIB) \
28 	$(FFTLIB) \
29 	$(PLUGINLIB) \
30 	$(SBLIB) \
31 	$(COLVARSLIB) \
32EOF_EOF
33patch -u -l -b -F 5 -N --suffix=.preplumed "./src/ComputeMgr.C" << \EOF_EOF
34--- ./src/ComputeMgr.C.preplumed
35+++ ./src/ComputeMgr.C
36@@ -92,12 +92,135 @@
37 #include "GlobalMasterSMD.h"
38 #include "GlobalMasterTMD.h"
39 #include "GlobalMasterSymmetry.h"
40 #include "GlobalMasterEasy.h"
41 #include "GlobalMasterMisc.h"
42+// PLUMED
43+#include "../Plumed.h"
44+// END PLUMED
45 #include "GlobalMasterFreeEnergy.h"
46 #include "GlobalMasterColvars.h"
47+class GlobalMasterPlumed:
48+  public  GlobalMasterEasy,
49+  private PLMD::Plumed
50+{
51+  std::vector<int> index;
52+  std::vector<double> positions;
53+  std::vector<double> forces;
54+  std::vector<double> masses;
55+  std::vector<double> charges;
56+  double box[3][3];
57+  SimParameters *spar;
58+public:
59+  GlobalMasterPlumed():
60+    GlobalMasterEasy("plumedScript")
61+  {
62+    easy_init(config);
63+  }
64+  void easy_init(const char* config)
65+  {
66+    int realPrecision=sizeof(double);
67+    cmd("setRealPrecision",&realPrecision);
68+    spar=Node::Object()->simParameters;
69+
70+    double energyUnits=4.184;
71+    double lengthUnits=0.1;
72+    double timeUnits=0.001;
73+    cmd("setMDEnergyUnits",&energyUnits);
74+    cmd("setMDLengthUnits",&lengthUnits);
75+    cmd("setMDTimeUnits",&timeUnits);
76+
77+    cmd("setPlumedDat",spar->plumedFilename);
78+
79+    int natoms=Node::Object()->molecule->numAtoms;
80+    cmd("setNatoms",&natoms);
81+
82+    cmd("setMDEngine","namd");
83+
84+    cmd("setLog",stdout);
85+
86+    double dt=spar->dt;
87+    cmd("setTimestep",&dt);
88+    cmd("setNoVirial");
89+    cmd("init");
90+
91+    int s=step+1;
92+    cmd("setStep",&s);
93+    share();
94+
95+  }
96+
97+  void share(){
98+    const int* p;
99+    int n;
100+    bool redo=false;
101+    cmd("prepareDependencies");
102+    cmd("createFullList",&n);
103+    cmd("getFullList",&p);
104+    if(index.size()!=n)redo=true;
105+    if(!redo) for(int i=0;i<n;i++) if(index[i]!=p[i]) { redo=true; break;};
106+    if(redo){
107+      index.resize(n);
108+      masses.resize(n);
109+      modifyRequestedAtoms().resize(0);
110+      for(int i=0;i<n;i++){
111+        requestAtom(p[i]);
112+        masses[i]=getMass(p[i]);
113+        index[i]=p[i];
114+      };
115+      positions.resize(3*n);
116+      forces.resize(3*n);
117+      charges.resize(n);
118+
119+      cmd("setAtomsNlocal",&n);
120+      cmd("setAtomsGatindex",&index[0]);
121+    }
122+    cmd("clearFullList");
123+  }
124+
125+  void easy_calc(){
126+
127+    for(int i=0;i<index.size();i++){
128+      Vector coord;
129+      getPosition(index[i],coord);
130+      positions[3*i+0]=coord.x;
131+      positions[3*i+1]=coord.y;
132+      positions[3*i+2]=coord.z;
133+      masses[i]=Node::Object()->molecule->atommass(index[i]);
134+      charges[i]=Node::Object()->molecule->atomcharge(index[i]);
135+    };
136+
137+    if(spar->lattice.volume()>0.0){
138+      for(int i=0;i<3;i++) box[0][i]=spar->lattice.a()[i];
139+      for(int i=0;i<3;i++) box[1][i]=spar->lattice.b()[i];
140+      for(int i=0;i<3;i++) box[2][i]=spar->lattice.c()[i];
141+      cmd("setBox",&box[0][0]);
142+    }
143+
144+    for(int i=0;i<forces.size();i++) forces[i]=0.0;
145+
146+    cmd("setMasses",&masses[0]); cmd("setCharges",&charges[0]);
147+    cmd("setPositions",&positions[0]);
148+    cmd("setForces",&forces[0]);
149+
150+    cmd("shareData");
151+    cmd("performCalc");
152+
153+    for(int i=0;i<index.size();i++){
154+      Vector f;
155+      f.x=forces[3*i+0];
156+      f.y=forces[3*i+1];
157+      f.z=forces[3*i+2];
158+      int j=addForce(index[i],f);
159+    };
160+
161+    int s=step+1;
162+    cmd("setStep",&s);
163+    share();
164+
165+  }
166+};
167
168 #include "ComputeNonbondedMICKernel.h"
169
170 #include "DeviceCUDA.h"
171 #ifdef NAMD_CUDA
172@@ -1048,10 +1171,15 @@
173             masterServerObject->addClient(new GlobalMasterSymmetry());
174         if (simParams->TMDOn)
175             masterServerObject->addClient(new GlobalMasterTMD());
176         if (simParams->miscForcesOn)
177             masterServerObject->addClient(new GlobalMasterMisc());
178+     // PLUMED
179+        if (simParams->plumedOn){
180+          masterServerObject->addClient(new GlobalMasterPlumed());
181+        }
182+     // END PLUMED
183         if ( simParams->freeEnergyOn )
184             masterServerObject->addClient(new GlobalMasterFreeEnergy());
185 		if ( simParams->colvarsOn )
186 			masterServerObject->addClient(new GlobalMasterColvars());
187
188EOF_EOF
189patch -u -l -b -F 5 -N --suffix=.preplumed "./src/SimParameters.C" << \EOF_EOF
190--- ./src/SimParameters.C.preplumed
191+++ ./src/SimParameters.C
192@@ -1778,10 +1778,17 @@
193    opts.optionalB("main", "GBIS", "Use GB implicit solvent?",
194       &GBISOn, FALSE);
195    opts.optionalB("main", "GBISSer", "Use GB implicit solvent?",
196       &GBISserOn, FALSE);
197
198+
199+     // PLUMED
200+     //// plumed
201+     opts.optionalB("main", "plumed", "Is PLUMED active?",&plumedOn,FALSE);
202+     opts.require("plumed", "plumedfile","PLUMED script",PARSE_STRING);
203+     // END PLUMED
204+
205    opts.optional("GBIS", "solventDielectric",
206       "Solvent Dielectric", &solvent_dielectric, 78.5);
207    opts.optional("GBIS", "intrinsicRadiusOffset",
208       "Coulomb Radius Offset", &coulomb_radius_offset, 0.09);
209    opts.optional("GBIS", "ionConcentration",
210@@ -5447,11 +5454,11 @@
211
212    // Global forces configuration
213
214    globalForcesOn = ( tclForcesOn || freeEnergyOn || miscForcesOn ||
215                       (IMDon && ! (IMDignore || IMDignoreForces)) || SMDOn || TMDOn ||
216-                      colvarsOn || symmetryOn || qmForcesOn );
217+                      colvarsOn || symmetryOn || qmForcesOn || plumedOn );
218
219
220    if (tclForcesOn)
221    {
222      iout << iINFO << "TCL GLOBAL FORCES ACTIVE\n";
223@@ -5532,10 +5539,25 @@
224      }
225
226      iout << endi;
227    }
228
229+   // PLUMED
230+   if (plumedOn)
231+   {
232+       iout << iINFO << "PLUMED ACTIVE\n";
233+
234+       current = config->find("plumedfile");
235+       iout << iINFO << "PLUMED CONFIG FILE   "<< current->data << "\n" << endi;
236+       strcpy(plumedFilename,current->data);
237+
238+       ifstream plumedFile(current->data);
239+       if ( ! plumedFile ) NAMD_die("Error reading PLUMED config file.\n");
240+
241+   }
242+   // END PLUMED
243+
244    if (IMDon)
245    {
246      iout << iINFO << "INTERACTIVE MD ACTIVE\n";
247      iout << iINFO << "INTERACTIVE MD PORT    " << IMDport << "\n";
248      iout << iINFO << "INTERACTIVE MD FREQ    " << IMDfreq << "\n";
249EOF_EOF
250patch -u -l -b -F 5 -N --suffix=.preplumed "./src/SimParameters.h" << \EOF_EOF
251--- ./src/SimParameters.h.preplumed
252+++ ./src/SimParameters.h
253@@ -547,10 +547,13 @@
254 	char tclBCArgs[128];		//  Extra args for calcforces command
255 	Bool freeEnergyOn;		//  Doing free energy perturbation?
256 	Bool miscForcesOn;		//  Using misc forces?
257 	Bool colvarsOn;         //  Using the colvars module?
258
259+        Bool plumedOn;                  //  Using Plumed?
260+        char plumedFilename[1024];      // Plumed filename
261+
262 	Bool fixedAtomsOn;		//  Are there fixed atoms?
263 	Bool fixedAtomsForces;		//  Calculate forces anyway?
264 	Bool fixedAtomsForceOutput; // Output fixed forces?
265
266 	Bool langevinOnAtStartup;	//  Ensure that langevin is set up properly
267EOF_EOF
268