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