Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2022.
3 :
4 : CVs originally developed by the Jochen Hub group from the University of Saarland (Germany)
5 : and adapted and implemented in PLUMED by Ary Lautaro Di Bartolo and Diego Masone from the
6 : National University of Cuyo (Argentina).
7 :
8 : The membranefusion module is free software: you can redistribute it and/or modify
9 : it under the terms of the GNU Lesser General Public License as published by
10 : the Free Software Foundation, either version 3 of the License, or
11 : (at your option) any later version.
12 :
13 : The membranefusion module is distributed in the hope that it will be useful,
14 : but WITHOUT ANY WARRANTY; without even the implied warranty of
15 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 : GNU Lesser General Public License for more details.
17 :
18 : You should have received a copy of the GNU Lesser General Public License
19 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
20 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
21 : #include "colvar/Colvar.h"
22 : #include "core/ActionRegister.h"
23 : #include <cmath>
24 : #ifdef _OPENMP
25 : #if _OPENMP >= 201307
26 : #include <omp.h>
27 : #endif
28 : #endif
29 :
30 : namespace PLMD {
31 : namespace membranefusion {
32 : //+PLUMEDOC MEMBRANEFUSIONMOD_COLVAR MEMFUSIONP
33 : /*
34 : Calculate a CV that can induce the formation of the hemifusion stalk between two initially flat and planar bilayers.
35 :
36 : Calculate the collective variable designed by Hub and collaborators \cite Hub2017 and
37 : implemented into PLUMED by Masone and collaborators \cite DiBartolo2022 .
38 : This CV is capable of inducing the formation of the hemifusion stalk between two initially flat and planar bilayers
39 : surrounded by water molecules.
40 :
41 : \f[
42 : \xi_f = \frac{1}{N_{sf}} \sum_{s=0}^{N_{sf}-1} \delta_{sf} (N_{sf}^{(p)})
43 : \f]
44 :
45 : Where \f$\xi_f\f$ is the CV, \f$N_{sf}\f$ is the number of slices of the cylinder that make up the CV,
46 : \f$\delta_{sf}\f$ is a continuos function in the interval [0 1] (\f$\delta_{sf} = 0\f$ for no beads in the slice s, and
47 : \f$\delta_{sf} = 1\f$ for 1 or more beads in the slice s) and \f$N_{sf}^{(p)}\f$ accounts for the number of tail beads
48 : within the slice s.
49 :
50 : \par Examples
51 :
52 : This example induces a hemifusion stalk (\f$\xi_f = 0.85\f$) from a pair of initially flat membranes (\f$\xi_f = 0.2\f$).
53 :
54 : \plumedfile
55 : lMem: GROUP ATOMS=1-12288 #All the lower membrane beads.
56 : uMem: GROUP ATOMS=12289-24576 #All the upper membrane beads.
57 : tails: GROUP ATOMS=8-24572:12,12-24576:12 #All the lipid tails beads (from the lower and upper membrane).
58 :
59 : memFusion: MEMFUSIONP UMEMBRANE=uMem LMEMBRANE=lMem TAILS=tails NSMEM=70 DSMEM=0.1 HMEM=0.25 RCYLMEM=1.75 ZETAMEM=0.5
60 :
61 : MOVINGRESTRAINT ...
62 : ARG=memFusion
63 : STEP0=0 AT0=0.2 KAPPA0=10000.0
64 : STEP1=500000 AT1=0.85 KAPPA1=10000.0
65 : ...
66 :
67 : PRINT ARG=memFusion FILE=COLVAR STRIDE=1
68 :
69 : \endplumedfile
70 :
71 : You can test this CV with another example in this <a href="https://github.com/lautarodibartolo/MemFusion/tree/main/ExampleParallel">GitHub folder</a>.
72 :
73 : */
74 : //+ENDPLUMEDOC
75 :
76 : class memFusionP : public Colvar {
77 : std::vector<AtomNumber> UMEM, LMEM, TAILS;
78 : std::vector<double> NSMEM, DSMEM, HMEM, RCYLMEM, ZETAMEM, ONEOVERS2C2CUTOFF, XCYL, YCYL;
79 :
80 : public:
81 : explicit memFusionP(const ActionOptions &);
82 : void calculate() override;
83 : static void registerKeywords(Keywords &keys);
84 : };
85 :
86 : PLUMED_REGISTER_ACTION(memFusionP, "MEMFUSIONP")
87 :
88 3 : void memFusionP::registerKeywords(Keywords &keys) {
89 3 : Colvar::registerKeywords(keys);
90 3 : keys.add("atoms", "UMEMBRANE", "all the beads of the upper membrane");
91 3 : keys.add("atoms", "LMEMBRANE", "all the beads of the lower membrane");
92 3 : keys.add("atoms", "TAILS", "all the tail beads of the system");
93 3 : keys.add("compulsory", "NSMEM", "the number of slices of the membrane fusion cylinder in such a way that when the bilayers are flat and parallel the CV is equal to 0.2.");
94 3 : keys.add("optional", "DSMEM", "( default=0.1) thickness of the slices of the membrane fusion cylinder.");
95 3 : keys.add("optional", "HMEM", "( default=0.25 ) parameter of the step function θ(x,h) for the membrane fusion.");
96 3 : keys.add("optional", "RCYLMEM", "( default=1.75 ) the radius of the membrane fusion cylinder.");
97 3 : keys.add("optional", "ZETAMEM", "( default=0.5 ) occupation factor.");
98 3 : keys.add("optional", "ONEOVERS2C2CUTOFF", "( default=500 ) cut off large values for the derivative of the atan2 function.");
99 3 : keys.add("optional", "XCYL", "X coordinate of the fixed cylinder, if not present this will be calculated.");
100 3 : keys.add("optional", "YCYL", "Y coordinate of the fixed cylinder, if not present this will be calculated.");
101 6 : keys.setValueDescription("scalar","the value of the CV");
102 3 : }
103 :
104 1 : memFusionP::memFusionP(const ActionOptions &ao) : PLUMED_COLVAR_INIT(ao) {
105 2 : parseAtomList("UMEMBRANE", UMEM);
106 1 : if (UMEM.size() == 0) {
107 0 : error("UMEMBRANE has not any atom specified.");
108 : }
109 :
110 2 : parseAtomList("LMEMBRANE", LMEM);
111 1 : if (LMEM.size() == 0) {
112 0 : error("LMEMBRANE has not any atom specified.");
113 : }
114 :
115 2 : parseAtomList("TAILS", TAILS);
116 1 : if (TAILS.size() == 0) {
117 0 : error("TAILS has not any atom specified.");
118 : }
119 :
120 2 : parseVector("NSMEM", NSMEM);
121 1 : if (NSMEM.size() > 1) {
122 0 : error("NSMEM cannot take more than one value.");
123 : }
124 :
125 2 : parseVector("DSMEM", DSMEM);
126 1 : if (DSMEM.size() > 1) {
127 0 : error("DSMEM cannot take more than one value.");
128 : }
129 1 : if (DSMEM.size() == 0) {
130 0 : DSMEM.push_back(0.1);
131 : }
132 :
133 2 : parseVector("HMEM", HMEM);
134 1 : if (HMEM.size() > 1) {
135 0 : error("HMEM cannot take more than one value.");
136 : }
137 1 : if (HMEM.size() == 0) {
138 0 : HMEM.push_back(0.25);
139 : }
140 :
141 2 : parseVector("RCYLMEM", RCYLMEM);
142 1 : if (RCYLMEM.size() > 1) {
143 0 : error("RCYLMEM cannot take more than one value.");
144 : }
145 1 : if (RCYLMEM.size() == 0) {
146 0 : RCYLMEM.push_back(1.75);
147 : }
148 :
149 2 : parseVector("ZETAMEM", ZETAMEM);
150 1 : if (ZETAMEM.size() > 1) {
151 0 : error("ZETA cannot take more than one value.");
152 : }
153 1 : if (ZETAMEM.size() == 0) {
154 0 : ZETAMEM.push_back(0.5);
155 : }
156 :
157 2 : parseVector("ONEOVERS2C2CUTOFF", ONEOVERS2C2CUTOFF);
158 1 : if (ONEOVERS2C2CUTOFF.size() > 1) {
159 0 : error("ONEOVERS2C2CUTOFF cannot take more than one value.");
160 : }
161 1 : if (ONEOVERS2C2CUTOFF.size() == 0) {
162 1 : ONEOVERS2C2CUTOFF.push_back(500);
163 : }
164 :
165 2 : parseVector("XCYL", XCYL);
166 1 : if (XCYL.size() > 1) {
167 0 : error("XCYL cannot take more than one value.");
168 : }
169 1 : if (XCYL.size() == 0) {
170 1 : XCYL.push_back(-1.0);
171 : }
172 :
173 2 : parseVector("YCYL", YCYL);
174 1 : if (YCYL.size() > 1) {
175 0 : error("YCYL cannot take more than one value.");
176 : }
177 1 : if (YCYL.size() == 0) {
178 1 : YCYL.push_back(-1.0);
179 : }
180 :
181 1 : checkRead();
182 :
183 : std::vector<AtomNumber> atoms;
184 12289 : for (unsigned i = 0; i < UMEM.size(); i++) {
185 12288 : atoms.push_back(UMEM[i]);
186 : }
187 12289 : for (unsigned i = 0; i < LMEM.size(); i++) {
188 12288 : atoms.push_back(LMEM[i]);
189 : }
190 4097 : for (unsigned i = 0; i < TAILS.size(); i++) {
191 4096 : atoms.push_back(TAILS[i]);
192 : }
193 :
194 1 : addValueWithDerivatives();
195 1 : setNotPeriodic();
196 1 : requestAtoms(atoms);
197 1 : }
198 :
199 3 : void memFusionP::calculate() {
200 : /**************************
201 : * *
202 : * System *
203 : * *
204 : **************************/
205 :
206 : // Box dimensions.
207 3 : double Lx = getBox()[0][0], Ly = getBox()[1][1], Lz = getBox()[2][2];
208 :
209 : // Z center of the upper membrane (uMem) and lower membrane (lMem) for systems with PBC: https://en.wikipedia.org/wiki/Center_of_mass#Systems_with_periodic_boundary_conditions .
210 : double ZuMem, ZuMemcos = 0.0, ZuMemsin = 0.0, uMemAngle, ZlMem, ZlMemcos = 0.0, ZlMemsin = 0.0, lMemAngle;
211 :
212 : #ifdef _OPENMP
213 : #if _OPENMP >= 201307
214 3 : #pragma omp parallel for private(uMemAngle, lMemAngle) reduction(+:ZuMemcos, ZuMemsin, ZlMemcos, ZlMemsin)
215 : #endif
216 : #endif
217 : for (unsigned i = 0; i < UMEM.size(); i++) {
218 : uMemAngle = 2.0 * M_PI * getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i)))[2];
219 : lMemAngle = 2.0 * M_PI * getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + UMEM.size())))[2];
220 : ZuMemcos += cos(uMemAngle);
221 : ZuMemsin += sin(uMemAngle);
222 : ZlMemcos += cos(lMemAngle);
223 : ZlMemsin += sin(lMemAngle);
224 : }
225 :
226 3 : ZuMemcos = ZuMemcos / UMEM.size();
227 3 : ZuMemsin = ZuMemsin / UMEM.size();
228 3 : ZlMemcos = ZlMemcos / UMEM.size();
229 3 : ZlMemsin = ZlMemsin / UMEM.size();
230 :
231 3 : ZuMem = Lz * (atan2(-ZuMemsin, -ZuMemcos) + M_PI) / (2.0 * M_PI);
232 3 : ZlMem = Lz * (atan2(-ZlMemsin, -ZlMemcos) + M_PI) / (2.0 * M_PI);
233 :
234 : // Z center of the boths membranes (upper and lower).
235 3 : double ZMems = (ZuMem + ZlMem) / 2.0;
236 :
237 : /*************************
238 : * *
239 : * Xi_Mem *
240 : * *
241 : **************************/
242 :
243 : // Quantity of beads of the membranes.
244 3 : unsigned membraneBeads = UMEM.size() + LMEM.size();
245 :
246 : // Z distance from the lipid tail to the geometric center of both membranes.
247 : double ZTailDistance;
248 :
249 : // Z position of the first slice.
250 3 : double firstSliceZ_Mem = ZMems + (0.0 + 0.5 - NSMEM[0] / 2.0) * DSMEM[0];
251 :
252 : // Z distance between the first slice and the Z center of the membrane.
253 6 : double firstSliceZDist_Mem = pbcDistance(Vector(0.0, 0.0, firstSliceZ_Mem), Vector(0.0, 0.0, ZMems))[2];
254 :
255 : // Position in the cylinder.
256 : double PositionS_Mem;
257 :
258 : // Slices to analyze per particle.
259 3 : std::vector<unsigned> s1_Mem(TAILS.size()), s2_Mem(TAILS.size());
260 :
261 : // Mark the particles to analyze.
262 3 : std::vector<double> analyzeThisParticle_Mem(TAILS.size());
263 :
264 : // Eq. 7 Hub & Awasthi JCTC 2017.
265 3 : std::vector<double> faxial_Mem(TAILS.size() * NSMEM[0]);
266 :
267 : // Eq. 16 Hub & Awasthi JCTC 2017.
268 3 : std::vector<double> d_faxial_Mem_dz(TAILS.size() * NSMEM[0]);
269 :
270 : // Eq. 10 Hub & Awasthi JCTC 2017.
271 3 : std::vector<double> Fs_Mem(NSMEM[0]);
272 :
273 : // Eq. 11 Hub & Awasthi JCTC 2017.
274 3 : std::vector<double> ws_Mem(NSMEM[0]);
275 :
276 : // Eq. 10 Hub & Awasthi JCTC 2017.
277 : double W_Mem = 0.0;
278 :
279 : // Eq. 21 and 22 Hub & Awasthi JCTC 2017.
280 3 : std::vector<double> sx_Mem(NSMEM[0]), sy_Mem(NSMEM[0]), cx_Mem(NSMEM[0]), cy_Mem(NSMEM[0]);
281 :
282 : // Eq. 10 Hub & Awasthi JCTC 2017.
283 : double Xsc_Mem = 0.0, Xcc_Mem = 0.0, Ysc_Mem = 0.0, Ycc_Mem = 0.0;
284 :
285 : // Aux.
286 : double x, aux;
287 :
288 : // Scaled position of the lipid tail respect the origin of coordinates.
289 3 : Vector TailPosition;
290 :
291 : // Thanks stack overflow.
292 : #ifdef _OPENMP
293 : #if _OPENMP >= 201307
294 : #pragma omp declare reduction(vec_double_plus : std::vector<double> : \
295 : std::transform(omp_out.begin(), omp_out.end(), omp_in.begin(), omp_out.begin(), std::plus<double>())) \
296 : initializer(omp_priv = decltype(omp_orig)(omp_orig.size()))
297 : #endif
298 : #endif
299 :
300 : #ifdef _OPENMP
301 : #if _OPENMP >= 201307
302 3 : #pragma omp parallel for private(ZTailDistance, PositionS_Mem, TailPosition, x, aux) reduction(vec_double_plus:Fs_Mem, sx_Mem, sy_Mem, cx_Mem, cy_Mem)
303 : #endif
304 : #endif
305 : for (unsigned i = 0; i < TAILS.size(); i++) {
306 : ZTailDistance = pbcDistance(Vector(0.0, 0.0, ZMems), getPosition(i + membraneBeads))[2];
307 : PositionS_Mem = (ZTailDistance + firstSliceZDist_Mem) / DSMEM[0];
308 : // If the following condition is met the particle is in the Z space of the cylinder.
309 : if ((PositionS_Mem >= (-0.5 - HMEM[0])) && (PositionS_Mem <= (NSMEM[0] + 0.5 - 1.0 + HMEM[0]))) {
310 : analyzeThisParticle_Mem[i] = 1.0;
311 : // Defining the slices to analyze each particle.
312 : if (PositionS_Mem < 1) {
313 : s1_Mem[i] = 0;
314 : s2_Mem[i] = 2;
315 : } else if (PositionS_Mem <= (NSMEM[0] - 2.0)) {
316 : s1_Mem[i] = floor(PositionS_Mem) - 1;
317 : s2_Mem[i] = floor(PositionS_Mem) + 1;
318 : } else {
319 : s1_Mem[i] = NSMEM[0] - 3;
320 : s2_Mem[i] = NSMEM[0] - 1;
321 : }
322 :
323 : TailPosition = getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + membraneBeads)));
324 :
325 : for (unsigned s = s1_Mem[i]; s <= s2_Mem[i]; s++) {
326 : x = (ZTailDistance - (s + 0.5 - NSMEM[0] / 2.0) * DSMEM[0]) * 2.0 / DSMEM[0];
327 : if (!((x <= -1.0 - HMEM[0]) || (x >= 1.0 + HMEM[0]))) {
328 : if (((-1.0 + HMEM[0]) <= x) && (x <= (1.0 - HMEM[0]))) {
329 : faxial_Mem[i + TAILS.size() * s] = 1.0;
330 : Fs_Mem[s] += 1.0;
331 : sx_Mem[s] += sin(2.0 * M_PI * TailPosition[0]);
332 : sy_Mem[s] += sin(2.0 * M_PI * TailPosition[1]);
333 : cx_Mem[s] += cos(2.0 * M_PI * TailPosition[0]);
334 : cy_Mem[s] += cos(2.0 * M_PI * TailPosition[1]);
335 : } else if (((1.0 - HMEM[0]) < x) && (x < (1.0 + HMEM[0]))) {
336 : aux = 0.5 - ((3.0 * x - 3.0) / (4.0 * HMEM[0])) + (pow((x - 1.0), 3) / (4.0 * pow(HMEM[0], 3)));
337 : faxial_Mem[i + TAILS.size() * s] = aux;
338 : d_faxial_Mem_dz[i + TAILS.size() * s] = ((-3.0 / (4.0 * HMEM[0])) + ((3.0 * pow((x - 1), 2)) / (4.0 * pow(HMEM[0], 3)))) * 2.0 / DSMEM[0];
339 : Fs_Mem[s] += aux;
340 : sx_Mem[s] += aux * sin(2.0 * M_PI * TailPosition[0]);
341 : sy_Mem[s] += aux * sin(2.0 * M_PI * TailPosition[1]);
342 : cx_Mem[s] += aux * cos(2.0 * M_PI * TailPosition[0]);
343 : cy_Mem[s] += aux * cos(2.0 * M_PI * TailPosition[1]);
344 : } else if (((-1.0 - HMEM[0]) < x) && (x < (-1.0 + HMEM[0]))) {
345 : aux = 0.5 + ((3.0 * x + 3.0) / (4.0 * HMEM[0])) - (pow((x + 1.0), 3) / (4.0 * pow(HMEM[0], 3)));
346 : faxial_Mem[i + TAILS.size() * s] = aux;
347 : d_faxial_Mem_dz[i + TAILS.size() * s] = ((3.0 / (4.0 * HMEM[0])) - ((3.0 * pow((x + 1), 2)) / (4.0 * pow(HMEM[0], 3)))) * 2.0 / DSMEM[0];
348 : Fs_Mem[s] += aux;
349 : sx_Mem[s] += (aux * sin(2.0 * M_PI * TailPosition[0]));
350 : sy_Mem[s] += (aux * sin(2.0 * M_PI * TailPosition[1]));
351 : cx_Mem[s] += (aux * cos(2.0 * M_PI * TailPosition[0]));
352 : cy_Mem[s] += (aux * cos(2.0 * M_PI * TailPosition[1]));
353 : }
354 : }
355 : }
356 : }
357 : }
358 :
359 213 : for (unsigned s = 0; s < NSMEM[0]; s++) {
360 210 : if (Fs_Mem[s] != 0.0) {
361 106 : ws_Mem[s] = tanh(Fs_Mem[s]);
362 106 : W_Mem += ws_Mem[s];
363 106 : sx_Mem[s] = sx_Mem[s] / Fs_Mem[s];
364 106 : sy_Mem[s] = sy_Mem[s] / Fs_Mem[s];
365 106 : cx_Mem[s] = cx_Mem[s] / Fs_Mem[s];
366 106 : cy_Mem[s] = cy_Mem[s] / Fs_Mem[s];
367 106 : Xsc_Mem += sx_Mem[s] * ws_Mem[s];
368 106 : Ysc_Mem += sy_Mem[s] * ws_Mem[s];
369 106 : Xcc_Mem += cx_Mem[s] * ws_Mem[s];
370 106 : Ycc_Mem += cy_Mem[s] * ws_Mem[s];
371 : }
372 : }
373 :
374 3 : Xsc_Mem = Xsc_Mem / W_Mem;
375 3 : Ysc_Mem = Ysc_Mem / W_Mem;
376 3 : Xcc_Mem = Xcc_Mem / W_Mem;
377 3 : Ycc_Mem = Ycc_Mem / W_Mem;
378 :
379 : // Eq. 12 Hub & Awasthi JCTC 2017.
380 : double Xcyl_Mem, Ycyl_Mem;
381 :
382 3 : if ((XCYL[0] > 0.0) && (YCYL[0] > 0.0)) {
383 : Xcyl_Mem = XCYL[0];
384 : Ycyl_Mem = YCYL[0];
385 : } else {
386 3 : Xcyl_Mem = (atan2(-Xsc_Mem, -Xcc_Mem) + M_PI) * Lx / (2 * M_PI);
387 3 : Ycyl_Mem = (atan2(-Ysc_Mem, -Ycc_Mem) + M_PI) * Ly / (2 * M_PI);
388 : }
389 :
390 : // Eq. 25, 26 and 27 Hub & Awasthi JCTC 2017.
391 : double d_sx_Mem_dx, d_sx_Mem_dz, d_sy_Mem_dy, d_sy_Mem_dz, d_cx_Mem_dx, d_cx_Mem_dz, d_cy_Mem_dy, d_cy_Mem_dz;
392 :
393 : // Eq. 29 Hub & Awasthi JCTC 2017.
394 : double d_ws_Mem_dz;
395 :
396 : // Eq. 31, 32 and 33 Hub & Awasthi JCTC 2017
397 : double d_Xsc_Mem_dx, d_Xsc_Mem_dz, d_Xcc_Mem_dx, d_Xcc_Mem_dz, d_Ysc_Mem_dy, d_Ysc_Mem_dz, d_Ycc_Mem_dy, d_Ycc_Mem_dz;
398 :
399 : // Center of the cylinder. XY components are calculated (or defined), Z is the Z geometric center of the membranes of the system.
400 6 : Vector xyzCyl_Mem = pbcDistance(Vector(0.0, 0.0, 0.0), Vector(Xcyl_Mem, Ycyl_Mem, ZMems));
401 :
402 : // Distances from the lipid tails to center of the cylinder.
403 3 : std::vector<Vector> CylDistances_Mem(TAILS.size());
404 :
405 : // XY distance from the lipid tails to the center of the cylinder.
406 : double ri_Mem;
407 :
408 : // Eq. 8 Hub & Awasthi JCTC 2017.
409 : double fradial_Mem = 0;
410 :
411 : // Eq. 15 Hub & Awasthi JCTC 2017.
412 3 : std::vector<double> d_fradial_Mem_dx(TAILS.size()), d_fradial_Mem_dy(TAILS.size());
413 :
414 : // Eq. 35, 36, 37 and 38 Hub & Awasthi JCTC 2017.
415 3 : std::vector<double> d_Xcyl_Mem_dx(TAILS.size()), d_Xcyl_Mem_dz(TAILS.size()), d_Ycyl_Mem_dy(TAILS.size()), d_Ycyl_Mem_dz(TAILS.size());
416 :
417 : // To avoid rare instabilities auxX_Mem and auxY_Mem are truncated at a configurable value (default = 500).
418 3 : double auxX_Mem = (1 / (pow(Xsc_Mem, 2) + pow(Xcc_Mem, 2))), auxY_Mem = (1 / (pow(Ysc_Mem, 2) + pow(Ycc_Mem, 2)));
419 :
420 3 : if (auxX_Mem > ONEOVERS2C2CUTOFF[0]) {
421 0 : auxX_Mem = Lx * ONEOVERS2C2CUTOFF[0] / (2 * M_PI);
422 : } else {
423 3 : auxX_Mem = Lx * auxX_Mem / (2 * M_PI);
424 : }
425 :
426 3 : if (auxY_Mem > ONEOVERS2C2CUTOFF[0]) {
427 0 : auxY_Mem = Ly * ONEOVERS2C2CUTOFF[0] / (2 * M_PI);
428 : } else {
429 3 : auxY_Mem = Ly * auxY_Mem / (2 * M_PI);
430 : }
431 :
432 : // Number of lipid tails within the slice s of the membranes cylinder.
433 3 : std::vector<double> Nsp_Mem(NSMEM[0]), psi_Mem(NSMEM[0]), d_psi_Mem(NSMEM[0]);
434 :
435 : // Eq. 3 Hub & Awasthi JCTC 2017.
436 3 : double b_Mem = (ZETAMEM[0] / (1.0 - ZETAMEM[0])), c_Mem = ((1.0 - ZETAMEM[0]) * exp(b_Mem));
437 :
438 : // Eq. 19 Hub & Awasthi JCTC 2017.
439 3 : std::vector<double> fradial_Mem_d_faxial_Mem_dz(TAILS.size() * NSMEM[0]);
440 :
441 : // Eq. 20 Hub & Awasthi JCTC 2017.
442 3 : std::vector<double> Axs_Mem(NSMEM[0]), Ays_Mem(NSMEM[0]);
443 :
444 : // Eq. 1 Hub & Awasthi JCTC 2017. This is the CV that describes de Pore Nucleation.
445 3 : double Xi_Mem = 0.0;
446 :
447 : #ifdef _OPENMP
448 : #if _OPENMP >= 201307
449 3 : #pragma omp parallel for private(TailPosition,d_Xsc_Mem_dx,d_Xcc_Mem_dx,d_Ysc_Mem_dy,d_Ycc_Mem_dy,d_Xsc_Mem_dz,d_Xcc_Mem_dz,d_Ysc_Mem_dz,d_Ycc_Mem_dz,d_sx_Mem_dx,d_sy_Mem_dy,d_cx_Mem_dx,d_cy_Mem_dy,d_sx_Mem_dz,d_sy_Mem_dz,d_cx_Mem_dz,d_cy_Mem_dz,d_ws_Mem_dz,ri_Mem,x,fradial_Mem) reduction(vec_double_plus: Nsp_Mem, Axs_Mem, Ays_Mem)
450 : #endif
451 : #endif
452 : for (unsigned i = 0; i < TAILS.size(); i++) {
453 : if (analyzeThisParticle_Mem[i]) {
454 : TailPosition = getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + membraneBeads)));
455 : d_Xsc_Mem_dx = 0.0;
456 : d_Xcc_Mem_dx = 0.0;
457 : d_Ysc_Mem_dy = 0.0;
458 : d_Ycc_Mem_dy = 0.0;
459 : d_Xsc_Mem_dz = 0.0;
460 : d_Xcc_Mem_dz = 0.0;
461 : d_Ysc_Mem_dz = 0.0;
462 : d_Ycc_Mem_dz = 0.0;
463 : for (unsigned s = s1_Mem[i]; s <= s2_Mem[i]; s++) {
464 : if (Fs_Mem[s] != 0.0) {
465 : d_sx_Mem_dx = faxial_Mem[i + TAILS.size() * s] * 2.0 * M_PI * cos(2.0 * M_PI * TailPosition[0]) / (Lx * Fs_Mem[s]);
466 : d_sy_Mem_dy = faxial_Mem[i + TAILS.size() * s] * 2.0 * M_PI * cos(2.0 * M_PI * TailPosition[1]) / (Ly * Fs_Mem[s]);
467 : d_cx_Mem_dx = -faxial_Mem[i + TAILS.size() * s] * 2.0 * M_PI * sin(2.0 * M_PI * TailPosition[0]) / (Lx * Fs_Mem[s]);
468 : d_cy_Mem_dy = -faxial_Mem[i + TAILS.size() * s] * 2.0 * M_PI * sin(2.0 * M_PI * TailPosition[1]) / (Ly * Fs_Mem[s]);
469 : d_Xsc_Mem_dx += ws_Mem[s] * d_sx_Mem_dx / W_Mem;
470 : d_Xcc_Mem_dx += ws_Mem[s] * d_cx_Mem_dx / W_Mem;
471 : d_Ysc_Mem_dy += ws_Mem[s] * d_sy_Mem_dy / W_Mem;
472 : d_Ycc_Mem_dy += ws_Mem[s] * d_cy_Mem_dy / W_Mem;
473 :
474 : d_sx_Mem_dz = d_faxial_Mem_dz[i + TAILS.size() * s] * (sin(2.0 * M_PI * TailPosition[0]) - sx_Mem[s]) / Fs_Mem[s];
475 : d_sy_Mem_dz = d_faxial_Mem_dz[i + TAILS.size() * s] * (sin(2.0 * M_PI * TailPosition[1]) - sy_Mem[s]) / Fs_Mem[s];
476 : d_cx_Mem_dz = d_faxial_Mem_dz[i + TAILS.size() * s] * (cos(2.0 * M_PI * TailPosition[0]) - cx_Mem[s]) / Fs_Mem[s];
477 : d_cy_Mem_dz = d_faxial_Mem_dz[i + TAILS.size() * s] * (cos(2.0 * M_PI * TailPosition[1]) - cy_Mem[s]) / Fs_Mem[s];
478 : d_ws_Mem_dz = (1 - pow(ws_Mem[s], 2)) * d_faxial_Mem_dz[i + TAILS.size() * s];
479 : d_Xsc_Mem_dz += (ws_Mem[s] * d_sx_Mem_dz + d_ws_Mem_dz * (sx_Mem[s] - Xsc_Mem)) / W_Mem;
480 : d_Xcc_Mem_dz += (ws_Mem[s] * d_cx_Mem_dz + d_ws_Mem_dz * (cx_Mem[s] - Xcc_Mem)) / W_Mem;
481 : d_Ysc_Mem_dz += (ws_Mem[s] * d_sy_Mem_dz + d_ws_Mem_dz * (sy_Mem[s] - Ysc_Mem)) / W_Mem;
482 : d_Ycc_Mem_dz += (ws_Mem[s] * d_cy_Mem_dz + d_ws_Mem_dz * (cy_Mem[s] - Ycc_Mem)) / W_Mem;
483 : }
484 : }
485 : d_Xcyl_Mem_dx[i] = auxX_Mem * (-Xsc_Mem * d_Xcc_Mem_dx + Xcc_Mem * d_Xsc_Mem_dx);
486 : d_Xcyl_Mem_dz[i] = auxX_Mem * (-Xsc_Mem * d_Xcc_Mem_dz + Xcc_Mem * d_Xsc_Mem_dz);
487 : d_Ycyl_Mem_dy[i] = auxY_Mem * (-Ysc_Mem * d_Ycc_Mem_dy + Ycc_Mem * d_Ysc_Mem_dy);
488 : d_Ycyl_Mem_dz[i] = auxY_Mem * (-Ysc_Mem * d_Ycc_Mem_dz + Ycc_Mem * d_Ysc_Mem_dz);
489 :
490 : CylDistances_Mem[i] = pbcDistance(xyzCyl_Mem, pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + membraneBeads)));
491 : ri_Mem = sqrt(pow(CylDistances_Mem[i][0], 2) + pow(CylDistances_Mem[i][1], 2));
492 : x = ri_Mem / RCYLMEM[0];
493 : if (!((x <= -1.0 - HMEM[0]) || (x >= 1.0 + HMEM[0]))) {
494 : if (((-1.0 + HMEM[0]) <= x) && (x <= (1.0 - HMEM[0]))) {
495 : fradial_Mem = 1.0;
496 : } else if (((1.0 - HMEM[0]) < x) && (x < (1.0 + HMEM[0]))) {
497 : fradial_Mem = 0.5 - ((3.0 * x - 3.0) / (4.0 * HMEM[0])) + (pow((x - 1.0), 3) / (4.0 * pow(HMEM[0], 3)));
498 : d_fradial_Mem_dx[i] = ((-3.0 / (4.0 * HMEM[0])) + ((3.0 * pow((x - 1), 2)) / (4.0 * pow(HMEM[0], 3)))) * CylDistances_Mem[i][0] / (RCYLMEM[0] * ri_Mem);
499 : d_fradial_Mem_dy[i] = ((-3.0 / (4.0 * HMEM[0])) + ((3.0 * pow((x - 1), 2)) / (4.0 * pow(HMEM[0], 3)))) * CylDistances_Mem[i][1] / (RCYLMEM[0] * ri_Mem);
500 : } else if (((-1.0 - HMEM[0]) < x) && (x < (-1.0 + HMEM[0]))) {
501 : fradial_Mem = 0.5 + ((3.0 * x + 3.0) / (4.0 * HMEM[0])) - (pow((x + 1.0), 3) / (4.0 * pow(HMEM[0], 3)));
502 : d_fradial_Mem_dx[i] = ((3.0 / (4.0 * HMEM[0])) - ((3.0 * pow((x + 1), 2)) / (4.0 * pow(HMEM[0], 3)))) * CylDistances_Mem[i][0] / (RCYLMEM[0] * ri_Mem);
503 : d_fradial_Mem_dy[i] = ((3.0 / (4.0 * HMEM[0])) - ((3.0 * pow((x + 1), 2)) / (4.0 * pow(HMEM[0], 3)))) * CylDistances_Mem[i][1] / (RCYLMEM[0] * ri_Mem);
504 : }
505 :
506 : for (unsigned s = s1_Mem[i]; s <= s2_Mem[i]; s++) {
507 : Nsp_Mem[s] += fradial_Mem * faxial_Mem[i + TAILS.size() * s];
508 : Axs_Mem[s] += faxial_Mem[i + TAILS.size() * s] * d_fradial_Mem_dx[i];
509 : Ays_Mem[s] += faxial_Mem[i + TAILS.size() * s] * d_fradial_Mem_dy[i];
510 : fradial_Mem_d_faxial_Mem_dz[i + TAILS.size() * s] = fradial_Mem * d_faxial_Mem_dz[i + TAILS.size() * s];
511 : }
512 : }
513 : }
514 : }
515 :
516 213 : for (unsigned s = 0; s < NSMEM[0]; s++) {
517 210 : if (Nsp_Mem[s] <= 1.0) {
518 166 : psi_Mem[s] = ZETAMEM[0] * Nsp_Mem[s];
519 166 : d_psi_Mem[s] = ZETAMEM[0];
520 166 : Xi_Mem += psi_Mem[s];
521 : } else {
522 44 : psi_Mem[s] = 1.0 - c_Mem * exp(-b_Mem * Nsp_Mem[s]);
523 44 : d_psi_Mem[s] = b_Mem * c_Mem * exp(-b_Mem * Nsp_Mem[s]);
524 44 : Xi_Mem += psi_Mem[s];
525 : }
526 : }
527 :
528 3 : Xi_Mem = Xi_Mem / NSMEM[0];
529 :
530 : // Eq. 18 Hub & Awasthi JCTC 2017.
531 3 : std::vector<double> faxial_Mem_d_fradial_Mem_dx(TAILS.size() * NSMEM[0]), faxial_Mem_d_fradial_Mem_dy(TAILS.size() * NSMEM[0]), faxial_Mem_d_fradial_Mem_dz(TAILS.size() * NSMEM[0]);
532 :
533 : // Eq. 13 Hub & Awasthi JCTC 2017.
534 3 : std::vector<Vector> derivatives_Mem(TAILS.size());
535 :
536 : #ifdef _OPENMP
537 : #if _OPENMP >= 201307
538 3 : #pragma omp parallel for private(aux)
539 : #endif
540 : #endif
541 : for (unsigned i = 0; i < TAILS.size(); i++) {
542 : if (analyzeThisParticle_Mem[i]) {
543 : for (unsigned s = s1_Mem[i]; s <= s2_Mem[i]; s++) {
544 : if (faxial_Mem[i + TAILS.size() * s]) {
545 : faxial_Mem_d_fradial_Mem_dx[i + TAILS.size() * s] = faxial_Mem[i + TAILS.size() * s] * d_fradial_Mem_dx[i] - d_Xcyl_Mem_dx[i] * Axs_Mem[s];
546 : faxial_Mem_d_fradial_Mem_dy[i + TAILS.size() * s] = faxial_Mem[i + TAILS.size() * s] * d_fradial_Mem_dy[i] - d_Ycyl_Mem_dy[i] * Ays_Mem[s];
547 : faxial_Mem_d_fradial_Mem_dz[i + TAILS.size() * s] = -d_Xcyl_Mem_dz[i] * Axs_Mem[s] - d_Ycyl_Mem_dz[i] * Ays_Mem[s];
548 : }
549 : }
550 :
551 : for (unsigned s = s1_Mem[i]; s <= s2_Mem[i]; s++) {
552 : aux = d_psi_Mem[s] / NSMEM[0];
553 : derivatives_Mem[i][0] += aux * faxial_Mem_d_fradial_Mem_dx[i + TAILS.size() * s];
554 : derivatives_Mem[i][1] += aux * faxial_Mem_d_fradial_Mem_dy[i + TAILS.size() * s];
555 : derivatives_Mem[i][2] += aux * (faxial_Mem_d_fradial_Mem_dz[i + TAILS.size() * s] + fradial_Mem_d_faxial_Mem_dz[i + TAILS.size() * s]);
556 : }
557 : }
558 : }
559 :
560 : // Derivatives and virial for the Xi_Mem.
561 3 : Tensor virial;
562 12291 : for (unsigned i = 0; i < TAILS.size(); i++) {
563 12288 : setAtomsDerivatives((i + membraneBeads), derivatives_Mem[i]);
564 12288 : virial -= Tensor(CylDistances_Mem[i], derivatives_Mem[i]);
565 : }
566 :
567 3 : setValue(Xi_Mem);
568 3 : setBoxDerivatives(virial);
569 3 : }
570 : }
571 : }
|
