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 : {
32 : namespace membranefusion
33 : {
34 : //+PLUMEDOC MEMBRANEFUSIONMOD_COLVAR FUSIONPORENUCLEATIONP
35 : /*
36 : A CV for inducing the nucleation of the fusion pore from a hemifusion stalk.
37 :
38 : Calculate the collective variable designed by Hub and collaborators \cite Hub2017 and
39 : implemented into PLUMED by Masone and collaborators.
40 : This CV is capable of inducing the nucleation of the fusion pore from a hemifusion stalk.
41 :
42 : \f[
43 : \xi_n = \frac{1}{N_{sn}} \sum_{s=0}^{N_{sn}-1} \delta_{sn} (N_{sn}^{(p)})
44 : \f]
45 :
46 : Where \f$\xi_n\f$ is the CV, \f$N_{sn}\f$ is the number of slices of the cylinder that make up the CV,
47 : \f$\delta_{sn}\f$ is a continuos function in the interval [0 1] (\f$\delta_{sf} = 0\f$ for no beads in the slice s, and
48 : \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 water and
49 : phosphateoxygens beads within the slice s.
50 :
51 : \par Examples
52 :
53 : This example induces the nucleation of the fusion pore (\f$\xi_n = 1.0\f$) from a hemifusion stalk (\f$\xi_n = 0.2\f$).
54 :
55 : \plumedfile
56 :
57 : lMem: GROUP ATOMS=1-10752,21505-22728,23953-24420 #All the lower membrane beads.
58 : uMem: GROUP ATOMS=10753-21504,22729-23952,24421-24888 #All the upper membrane beads.
59 : tails: GROUP ATOMS=8-23948:12,12-23952:12,23966-24884:18,23970-24888:18 #All the lipid tails beads (from the lower and upper membrane).
60 : waters: GROUP ATOMS=24889-56490 #All the water beads.
61 : po4: GROUP ATOMS=2-23942:12,23957-24875:18 #All the lipid phosphateoxygens beads.
62 :
63 : fusionPoreNucleation: FUSIONPORENUCLEATIONP UMEMBRANE=uMem LMEMBRANE=lMem TAILS=tails WATERS=waters PHOSPHATEOXYGENS=po4 NSMEM=85 NS=45
64 :
65 : MOVINGRESTRAINT ...
66 : ARG=fusionPoreNucleation
67 : STEP0=0 AT0=0.2 KAPPA0=10000.0
68 : STEP1=500000 AT1=1.0 KAPPA1=10000.0
69 : ...
70 :
71 : PRINT ARG=fusionPoreNucleation FILE=COLVAR STRIDE=1
72 :
73 : \endplumedfile
74 :
75 : */
76 : //+ENDPLUMEDOC
77 : class fusionPoreNucleationP : public Colvar
78 : {
79 : std::vector<AtomNumber> UMEM, LMEM, TAILS, WATERS, POXYGENS;
80 : std::vector<double> NSMEM, DSMEM, HMEM, NS, DS, HCH, RCYL, ZETA, ONEOVERS2C2CUTOFF, XCYL, YCYL;
81 :
82 : public:
83 : explicit fusionPoreNucleationP(const ActionOptions &);
84 : void calculate() override;
85 : static void registerKeywords(Keywords &keys);
86 : };
87 :
88 : PLUMED_REGISTER_ACTION(fusionPoreNucleationP, "FUSIONPORENUCLEATIONP")
89 :
90 3 : void fusionPoreNucleationP::registerKeywords(Keywords &keys)
91 : {
92 3 : Colvar::registerKeywords(keys);
93 6 : keys.add("atoms", "UMEMBRANE", "all the beads of the upper membrane.");
94 6 : keys.add("atoms", "LMEMBRANE", "all the beads of the lower membrane.");
95 6 : keys.add("atoms", "TAILS", "all the tail beads of the system.");
96 6 : keys.add("atoms", "WATERS", "all the water beads of the system.");
97 6 : keys.add("atoms", "PHOSPHATEOXYGENS", "all the lipid phosphateoxygens beads of the system.");
98 6 : keys.add("compulsory", "NSMEM", "the number of slices of the membrane fusion cylinder.");
99 6 : keys.add("optional", "DSMEM", "( default=0.1 ) thickness of the slices of the membrane fusion cylinder.");
100 6 : keys.add("optional", "HMEM", "( default=0.25 ) parameter of the step function θ(x,h) for the membrane fusion.");
101 6 : keys.add("compulsory", "NS", "the number of slices of the membrane-spanning cylinder in such a way that when the bilayers are flat and parallel the CV is equal to 0.2.");
102 6 : keys.add("optional", "DS", "( default=0.25 ) thickness of the slices of the membrane-spanning cylinder.");
103 6 : keys.add("optional", "HCH", "( default=0.25 ) parameter of the step function θ(x,h) for the CV.");
104 6 : keys.add("optional", "RCYL", "( default=0.8 ) the radius of the membrane-spanning cylinder.");
105 6 : keys.add("optional", "ZETA", "( default=0.75 ) parameter of the switch function ψ(x,ζ).");
106 6 : keys.add("optional", "ONEOVERS2C2CUTOFF", "( default=500 ) cut off large values for the derivative of the atan2 function to avoid violate energy.");
107 6 : keys.add("optional", "XCYL", "X coordinate of the fixed cylinder, if not present this will be calculated.");
108 6 : keys.add("optional", "YCYL", "X coordinate of the fixed cylinder, if not present this will be calculated.");
109 3 : keys.setValueDescription("the value of the CV");
110 3 : }
111 :
112 1 : fusionPoreNucleationP::fusionPoreNucleationP(const ActionOptions &ao) : PLUMED_COLVAR_INIT(ao)
113 : {
114 2 : parseAtomList("UMEMBRANE", UMEM);
115 1 : if (UMEM.size() == 0)
116 0 : error("UMEMBRANE has not any atom specified.");
117 :
118 2 : parseAtomList("LMEMBRANE", LMEM);
119 1 : if (LMEM.size() == 0)
120 0 : error("LMEMBRANE has not any atom specified.");
121 :
122 2 : parseAtomList("TAILS", TAILS);
123 1 : if (TAILS.size() == 0)
124 0 : error("TAILS has not any atom specified.");
125 :
126 2 : parseAtomList("WATERS", WATERS);
127 1 : if (WATERS.size() == 0)
128 0 : error("WATERS has not any atom specified.");
129 :
130 2 : parseAtomList("PHOSPHATEOXYGENS", POXYGENS);
131 1 : if (POXYGENS.size() == 0)
132 0 : error("PHOSPHATEOXYGENS has not any atom specified.");
133 :
134 2 : parseVector("NSMEM", NSMEM);
135 1 : if (NSMEM.size() > 1)
136 0 : error("NSMEM cannot take more than one value.");
137 :
138 2 : parseVector("DSMEM", DSMEM);
139 1 : if (DSMEM.size() > 1)
140 0 : error("DSMEM cannot take more than one value.");
141 1 : if (DSMEM.size() == 0)
142 0 : DSMEM.push_back(0.1);
143 :
144 2 : parseVector("HMEM", HMEM);
145 1 : if (HMEM.size() > 1)
146 0 : error("HMEM cannot take more than one value.");
147 1 : if (HMEM.size() == 0)
148 0 : HMEM.push_back(0.25);
149 :
150 2 : parseVector("NS", NS);
151 1 : if (NS.size() > 1)
152 0 : error("NS cannot take more than one value.");
153 :
154 2 : parseVector("DS", DS);
155 1 : if (DS.size() > 1)
156 0 : error("DS cannot take more than one value.");
157 1 : if (DS.size() == 0)
158 0 : DS.push_back(0.25);
159 :
160 2 : parseVector("HCH", HCH);
161 1 : if (HCH.size() > 1)
162 0 : error("H cannot take more than one value.");
163 1 : if (HCH.size() == 0)
164 0 : HCH.push_back(0.25);
165 :
166 2 : parseVector("RCYL", RCYL);
167 1 : if (RCYL.size() > 1)
168 0 : error("RCYL cannot take more than one value.");
169 1 : if (RCYL.size() == 0)
170 0 : RCYL.push_back(0.8);
171 :
172 2 : parseVector("ZETA", ZETA);
173 1 : if (ZETA.size() > 1)
174 0 : error("ZETA cannot take more than one value.");
175 1 : if (ZETA.size() == 0)
176 0 : ZETA.push_back(0.75);
177 :
178 2 : parseVector("ONEOVERS2C2CUTOFF", ONEOVERS2C2CUTOFF);
179 1 : if (ONEOVERS2C2CUTOFF.size() > 1)
180 0 : error("ONEOVERS2C2CUTOFF cannot take more than one value.");
181 1 : if (ONEOVERS2C2CUTOFF.size() == 0)
182 1 : ONEOVERS2C2CUTOFF.push_back(500);
183 :
184 2 : parseVector("XCYL", XCYL);
185 1 : if (XCYL.size() > 1)
186 0 : error("XCYL cannot take more than one value.");
187 1 : if (XCYL.size() == 0)
188 1 : XCYL.push_back(-1.0);
189 :
190 2 : parseVector("YCYL", YCYL);
191 1 : if (YCYL.size() > 1)
192 0 : error("YCYL cannot take more than one value.");
193 1 : if (YCYL.size() == 0)
194 1 : YCYL.push_back(-1.0);
195 :
196 1 : checkRead();
197 :
198 : std::vector<AtomNumber> atoms;
199 12445 : for (unsigned i = 0; i < UMEM.size(); i++)
200 : {
201 12444 : atoms.push_back(UMEM[i]);
202 : }
203 12445 : for (unsigned i = 0; i < LMEM.size(); i++)
204 : {
205 12444 : atoms.push_back(LMEM[i]);
206 : }
207 4097 : for (unsigned i = 0; i < TAILS.size(); i++)
208 : {
209 4096 : atoms.push_back(TAILS[i]);
210 : }
211 31603 : for (unsigned i = 0; i < WATERS.size(); i++)
212 : {
213 31602 : atoms.push_back(WATERS[i]);
214 : }
215 2049 : for (unsigned i = 0; i < POXYGENS.size(); i++)
216 : {
217 2048 : atoms.push_back(POXYGENS[i]);
218 : }
219 :
220 1 : addValueWithDerivatives();
221 1 : setNotPeriodic();
222 1 : requestAtoms(atoms);
223 1 : }
224 :
225 4 : void fusionPoreNucleationP::calculate()
226 : {
227 : /*************************
228 : * *
229 : * System *
230 : * *
231 : **************************/
232 :
233 : // Box dimensions.
234 4 : double Lx = getBox()[0][0], Ly = getBox()[1][1], Lz = getBox()[2][2];
235 :
236 : // 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 .
237 : double ZuMem, ZuMemcos = 0.0, ZuMemsin = 0.0, uMemAngle, ZlMem, ZlMemcos = 0.0, ZlMemsin = 0.0, lMemAngle;
238 :
239 : #ifdef _OPENMP
240 : #if _OPENMP >= 201307
241 4 : #pragma omp parallel for private(uMemAngle, lMemAngle) reduction(+:ZuMemcos, ZuMemsin, ZlMemcos, ZlMemsin)
242 : #endif
243 : #endif
244 : for (unsigned i = 0; i < UMEM.size(); i++)
245 : {
246 : uMemAngle = 2.0 * M_PI * getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i)))[2];
247 : lMemAngle = 2.0 * M_PI * getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + UMEM.size())))[2];
248 : ZuMemcos += cos(uMemAngle);
249 : ZuMemsin += sin(uMemAngle);
250 : ZlMemcos += cos(lMemAngle);
251 : ZlMemsin += sin(lMemAngle);
252 : }
253 4 : ZuMemcos = ZuMemcos / UMEM.size();
254 4 : ZuMemsin = ZuMemsin / UMEM.size();
255 4 : ZuMem = Lz * (atan2(-ZuMemsin, -ZuMemcos) + M_PI) / (2.0 * M_PI);
256 4 : ZlMemcos = ZlMemcos / UMEM.size();
257 4 : ZlMemsin = ZlMemsin / UMEM.size();
258 4 : ZlMem = Lz * (atan2(-ZlMemsin, -ZlMemcos) + M_PI) / (2.0 * M_PI);
259 :
260 : // Z center of the boths membranes (upper and lower).
261 4 : double ZMems = (ZuMem + ZlMem) / 2.0;
262 :
263 : /**************************
264 : * *
265 : * Xcyl_Mem & Ycyl_Mem *
266 : * *
267 : ***************************/
268 :
269 : // Quantity of beads of the membranes.
270 4 : unsigned membraneBeads = UMEM.size() + LMEM.size();
271 :
272 : // Z distance from the lipid tail to the geometric center of both membranes.
273 : double ZTailDistance;
274 :
275 : // Z position of the first slice.
276 4 : double firstSliceZ_Mem = ZMems + (0.0 + 0.5 - NSMEM[0] / 2.0) * DSMEM[0];
277 :
278 : // Z distance between the first slice and the Z center of the membrane.
279 8 : double firstSliceZDist_Mem = pbcDistance(Vector(0.0, 0.0, firstSliceZ_Mem), Vector(0.0, 0.0, ZMems))[2];
280 :
281 : // Position in the cylinder.
282 : double PositionS_Mem;
283 :
284 : // Slices to analyze per particle.
285 : unsigned s1_Mem, s2_Mem;
286 :
287 : // Eq. 7 Hub & Awasthi JCTC 2017.
288 4 : std::vector<double> faxial_Mem(TAILS.size() * NSMEM[0]);
289 :
290 : // Eq. 10 Hub & Awasthi JCTC 2017.
291 4 : std::vector<double> Fs_Mem(NSMEM[0]);
292 :
293 : // Eq. 11 Hub & Awasthi JCTC 2017.
294 4 : std::vector<double> ws_Mem(NSMEM[0]);
295 :
296 : // Eq. 10 Hub & Awasthi JCTC 2017.
297 : double W_Mem = 0.0;
298 :
299 : // Eq. 21 and 22 Hub & Awasthi JCTC 2017.
300 4 : std::vector<double> sx_Mem(NSMEM[0]), sy_Mem(NSMEM[0]), cx_Mem(NSMEM[0]), cy_Mem(NSMEM[0]);
301 :
302 : // Eq. 10 Hub & Awasthi JCTC 2017.
303 : double Xsc_Mem = 0.0, Xcc_Mem = 0.0, Ysc_Mem = 0.0, Ycc_Mem = 0.0;
304 :
305 : // Aux.
306 : double x, aux;
307 :
308 : // Scaled position of the lipid tail respect the origin of coordinates.
309 4 : Vector TailPosition;
310 :
311 : #ifdef _OPENMP
312 : #if _OPENMP >= 201307
313 : #pragma omp declare reduction(vec_double_plus : std::vector<double> : \
314 : std::transform(omp_out.begin(), omp_out.end(), omp_in.begin(), omp_out.begin(), std::plus<double>())) \
315 : initializer(omp_priv = decltype(omp_orig)(omp_orig.size()))
316 : #endif
317 : #endif
318 :
319 : #ifdef _OPENMP
320 : #if _OPENMP >= 201307
321 4 : #pragma omp parallel for private(ZTailDistance, PositionS_Mem, s1_Mem, s2_Mem, TailPosition, x, aux) reduction(vec_double_plus:Fs_Mem, sx_Mem, sy_Mem, cx_Mem, cy_Mem)
322 : #endif
323 : #endif
324 : for (unsigned i = 0; i < TAILS.size(); i++)
325 : {
326 : ZTailDistance = pbcDistance(Vector(0.0, 0.0, ZMems), getPosition(i + membraneBeads))[2];
327 : PositionS_Mem = (ZTailDistance + firstSliceZDist_Mem) / DSMEM[0];
328 : // If the following condition is met the particle is in the Z space of the cylinder.
329 : if ((PositionS_Mem >= (-0.5 - HMEM[0])) && (PositionS_Mem <= (NSMEM[0] + 0.5 - 1.0 + HMEM[0])))
330 : {
331 : //Defining the slices to analyze each particle.
332 : if (PositionS_Mem < 1)
333 : {
334 : s1_Mem = 0;
335 : s2_Mem = 2;
336 : }
337 : else if (PositionS_Mem <= (NSMEM[0] - 2.0))
338 : {
339 : s1_Mem = floor(PositionS_Mem) - 1;
340 : s2_Mem = floor(PositionS_Mem) + 1;
341 : }
342 : else
343 : {
344 : s1_Mem = NSMEM[0] - 3;
345 : s2_Mem = NSMEM[0] - 1;
346 : }
347 :
348 : TailPosition = getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + membraneBeads)));
349 :
350 : for (unsigned s = s1_Mem; s <= s2_Mem; s++)
351 : {
352 : x = (ZTailDistance - (s + 0.5 - NSMEM[0] / 2.0) * DSMEM[0]) * 2.0 / DSMEM[0];
353 : if (!((x <= -1.0 - HMEM[0]) || (x >= 1.0 + HMEM[0])))
354 : {
355 : if (((-1.0 + HMEM[0]) <= x) && (x <= (1.0 - HMEM[0])))
356 : {
357 : faxial_Mem[i + TAILS.size() * s] = 1.0;
358 : Fs_Mem[s] += 1.0;
359 : sx_Mem[s] += sin(2.0 * M_PI * TailPosition[0]);
360 : sy_Mem[s] += sin(2.0 * M_PI * TailPosition[1]);
361 : cx_Mem[s] += cos(2.0 * M_PI * TailPosition[0]);
362 : cy_Mem[s] += cos(2.0 * M_PI * TailPosition[1]);
363 : }
364 : else if (((1.0 - HMEM[0]) < x) && (x < (1.0 + HMEM[0])))
365 : {
366 : aux = 0.5 - ((3.0 * x - 3.0) / (4.0 * HMEM[0])) + (pow((x - 1.0), 3) / (4.0 * pow(HMEM[0], 3)));
367 : faxial_Mem[i + TAILS.size() * s] = aux;
368 : Fs_Mem[s] += aux;
369 : sx_Mem[s] += aux * sin(2.0 * M_PI * TailPosition[0]);
370 : sy_Mem[s] += aux * sin(2.0 * M_PI * TailPosition[1]);
371 : cx_Mem[s] += aux * cos(2.0 * M_PI * TailPosition[0]);
372 : cy_Mem[s] += aux * cos(2.0 * M_PI * TailPosition[1]);
373 : }
374 : else if (((-1.0 - HMEM[0]) < x) && (x < (-1.0 + HMEM[0])))
375 : {
376 : aux = 0.5 + ((3.0 * x + 3.0) / (4.0 * HMEM[0])) - (pow((x + 1.0), 3) / (4.0 * pow(HMEM[0], 3)));
377 : faxial_Mem[i + TAILS.size() * s] = aux;
378 : Fs_Mem[s] += aux;
379 : sx_Mem[s] += (aux * sin(2.0 * M_PI * TailPosition[0]));
380 : sy_Mem[s] += (aux * sin(2.0 * M_PI * TailPosition[1]));
381 : cx_Mem[s] += (aux * cos(2.0 * M_PI * TailPosition[0]));
382 : cy_Mem[s] += (aux * cos(2.0 * M_PI * TailPosition[1]));
383 : }
384 : }
385 : }
386 : }
387 : }
388 :
389 344 : for (unsigned s = 0; s < NSMEM[0]; s++)
390 : {
391 340 : if (Fs_Mem[s] != 0.0)
392 : {
393 339 : ws_Mem[s] = tanh(Fs_Mem[s]);
394 339 : W_Mem += ws_Mem[s];
395 339 : sx_Mem[s] = sx_Mem[s] / Fs_Mem[s];
396 339 : sy_Mem[s] = sy_Mem[s] / Fs_Mem[s];
397 339 : cx_Mem[s] = cx_Mem[s] / Fs_Mem[s];
398 339 : cy_Mem[s] = cy_Mem[s] / Fs_Mem[s];
399 339 : Xsc_Mem += sx_Mem[s] * ws_Mem[s];
400 339 : Ysc_Mem += sy_Mem[s] * ws_Mem[s];
401 339 : Xcc_Mem += cx_Mem[s] * ws_Mem[s];
402 339 : Ycc_Mem += cy_Mem[s] * ws_Mem[s];
403 : }
404 : }
405 :
406 4 : Xsc_Mem = Xsc_Mem / W_Mem;
407 4 : Ysc_Mem = Ysc_Mem / W_Mem;
408 4 : Xcc_Mem = Xcc_Mem / W_Mem;
409 4 : Ycc_Mem = Ycc_Mem / W_Mem;
410 :
411 : // Eq. 12 Hub & Awasthi JCTC 2017.
412 : double Xcyl_Mem, Ycyl_Mem;
413 :
414 4 : if ((XCYL[0] > 0.0) && (YCYL[0] > 0.0))
415 : {
416 : Xcyl_Mem = XCYL[0];
417 : Ycyl_Mem = YCYL[0];
418 : }
419 : else
420 : {
421 4 : Xcyl_Mem = (atan2(-Xsc_Mem, -Xcc_Mem) + M_PI) * Lx / (2 * M_PI);
422 4 : Ycyl_Mem = (atan2(-Ysc_Mem, -Ycc_Mem) + M_PI) * Ly / (2 * M_PI);
423 : }
424 :
425 : /*************************
426 : * *
427 : * Xi_n *
428 : * *
429 : **************************/
430 :
431 : // Eq. 1 Hub & Awasthi JCTC 2017. This is the CV that describes de Pore Nucleation.
432 4 : double Xi_n = 0.0;
433 :
434 : // Quantity of beads that could participate in the calculation of the Xi_Chain
435 4 : unsigned chainBeads = WATERS.size() + POXYGENS.size();
436 :
437 : // Quantity of beads that don't participate in the calculation of the Xi_Chain
438 4 : unsigned noChainBeads = (UMEM.size() * 2) + TAILS.size();
439 :
440 : // Z Distances from the oxygens to the geometric center of the membranes.
441 : double ZMemDistance;
442 :
443 : // Scaled positions of the oxygens to respect of the origin of coordinates.
444 4 : Vector Position;
445 :
446 : // Distance from the water/phosphate group to the defect cylinder.
447 4 : Vector distCylinder;
448 :
449 : // Center of the cylinder. XY components are calculated (or defined), Z is the Z geometric center of the membranes of the system.
450 8 : Vector xyzCyl_Mem = pbcDistance(Vector(0.0, 0.0, 0.0), Vector(Xcyl_Mem, Ycyl_Mem, ZMems));
451 :
452 : // Average of the radius of the water and lipid cylinder.
453 4 : double RCYLAVERAGE = RCYL[0] * (1 + HCH[0]);
454 :
455 : // Conditions.
456 : bool condition1, condition2, condition3;
457 :
458 : // Z position of the first slice.
459 4 : double firstSliceZ = ZMems + (0.0 + 0.5 - NS[0] / 2.0) * DS[0];
460 :
461 : // Z distance between the first slice and the Z center of the membrane.
462 4 : double firstSliceZDist = pbcDistance(Vector(0.0, 0.0, firstSliceZ), Vector(0.0, 0.0, ZMems))[2];
463 :
464 : // Position in the cylinder.
465 : double PositionS;
466 :
467 : // Mark the particles to analyze.
468 4 : std::vector<double> analyzeThisParticle(chainBeads);
469 :
470 : // Slices to analyze per particle.
471 4 : std::vector<unsigned> s1(chainBeads), s2(chainBeads);
472 :
473 : // Eq. 7 Hub & Awasthi JCTC 2017.
474 4 : std::vector<double> faxial(chainBeads * NS[0]);
475 :
476 : // Eq. 16 Hub & Awasthi JCTC 2017.
477 4 : std::vector<double> d_faxial_dz(chainBeads * NS[0]);
478 :
479 : // Eq. 10 Hub & Awasthi JCTC 2017.
480 4 : std::vector<double> Fs(NS[0]);
481 :
482 : // Eq. 11 Hub & Awasthi JCTC 2017.
483 4 : std::vector<double> ws(NS[0]);
484 :
485 : // Eq. 10 Hub & Awasthi JCTC 2017.
486 : double W = 0.0;
487 :
488 : // Eq. 21 and 22 Hub & Awasthi JCTC 2017.
489 4 : std::vector<double> sx(NS[0]), sy(NS[0]), cx(NS[0]), cy(NS[0]);
490 :
491 : // Eq. 10 Hub & Awasthi JCTC 2017.
492 : double Xsc = 0.0, Xcc = 0.0, Ysc = 0.0, Ycc = 0.0;
493 :
494 : #ifdef _OPENMP
495 : #if _OPENMP >= 201307
496 4 : #pragma omp parallel for private(distCylinder, aux, condition1, condition2, condition3, ZMemDistance, PositionS, Position, x) reduction(vec_double_plus:Fs, sx, sy, cx, cy)
497 : #endif
498 : #endif
499 : for (unsigned i = 0; i < chainBeads; i++)
500 : {
501 : distCylinder = pbcDistance(xyzCyl_Mem, pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + noChainBeads)));
502 : aux = sqrt(pow(distCylinder[0], 2) + pow(distCylinder[1], 2));
503 : condition1 = ((aux / RCYLAVERAGE) < 1.0);
504 : condition2 = ((pbcDistance(Vector(0.0, 0.0, ZuMem), getPosition(i + noChainBeads))[2] > 0) && (aux / RCYLAVERAGE) < 2.0);
505 : condition3 = ((pbcDistance(getPosition(i + noChainBeads), Vector(0.0, 0.0, ZlMem))[2] > 0) && (aux / RCYLAVERAGE) < 2.0);
506 : if (condition1 || condition2 || condition3)
507 : {
508 : ZMemDistance = pbcDistance(Vector(0.0, 0.0, ZMems), getPosition(i + noChainBeads))[2];
509 : PositionS = (ZMemDistance + firstSliceZDist) / DS[0];
510 : // If the following condition is met the particle is in the Z space of the cylinder.
511 : if ((PositionS >= (-0.5 - HCH[0])) && (PositionS <= (NS[0] + 0.5 - 1.0 + HCH[0])))
512 : {
513 : analyzeThisParticle[i] = 1.0;
514 :
515 : //Defining the slices to analyze each particle.
516 : if (PositionS < 1)
517 : {
518 : s1[i] = 0;
519 : s2[i] = 2;
520 : }
521 : else if (PositionS <= (NS[0] - 2.0))
522 : {
523 : s1[i] = floor(PositionS) - 1;
524 : s2[i] = floor(PositionS) + 1;
525 : }
526 : else
527 : {
528 : s1[i] = NS[0] - 3;
529 : s2[i] = NS[0] - 1;
530 : }
531 :
532 : Position = getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + noChainBeads)));
533 :
534 : for (unsigned s = s1[i]; s <= s2[i]; s++)
535 : {
536 : x = (ZMemDistance - (s + 0.5 - NS[0] / 2.0) * DS[0]) * 2.0 / DS[0];
537 : if (!((x <= -1.0 - HCH[0]) || (x >= 1.0 + HCH[0])))
538 : {
539 : if (((-1.0 + HCH[0]) <= x) && (x <= (1.0 - HCH[0])))
540 : {
541 : faxial[i + chainBeads * s] = 1.0;
542 : Fs[s] += 1.0;
543 : sx[s] += sin(2.0 * M_PI * Position[0]);
544 : sy[s] += sin(2.0 * M_PI * Position[1]);
545 : cx[s] += cos(2.0 * M_PI * Position[0]);
546 : cy[s] += cos(2.0 * M_PI * Position[1]);
547 : }
548 : else if (((1.0 - HCH[0]) < x) && (x < (1.0 + HCH[0])))
549 : {
550 : aux = 0.5 - ((3.0 * x - 3.0) / (4.0 * HCH[0])) + (pow((x - 1.0), 3) / (4.0 * pow(HCH[0], 3)));
551 : faxial[i + chainBeads * s] = aux;
552 : d_faxial_dz[i + chainBeads * s] = ((-3.0 / (4.0 * HCH[0])) + ((3.0 * pow((x - 1), 2)) / (4.0 * pow(HCH[0], 3)))) * 2.0 / DS[0];
553 : Fs[s] += aux;
554 : sx[s] += aux * sin(2.0 * M_PI * Position[0]);
555 : sy[s] += aux * sin(2.0 * M_PI * Position[1]);
556 : cx[s] += aux * cos(2.0 * M_PI * Position[0]);
557 : cy[s] += aux * cos(2.0 * M_PI * Position[1]);
558 : }
559 : else if (((-1.0 - HCH[0]) < x) && (x < (-1.0 + HCH[0])))
560 : {
561 : aux = 0.5 + ((3.0 * x + 3.0) / (4.0 * HCH[0])) - (pow((x + 1.0), 3) / (4.0 * pow(HCH[0], 3)));
562 : faxial[i + chainBeads * s] = aux;
563 : d_faxial_dz[i + chainBeads * s] = ((3.0 / (4.0 * HCH[0])) - ((3.0 * pow((x + 1), 2)) / (4.0 * pow(HCH[0], 3)))) * 2.0 / DS[0];
564 : Fs[s] += aux;
565 : sx[s] += (aux * sin(2.0 * M_PI * Position[0]));
566 : sy[s] += (aux * sin(2.0 * M_PI * Position[1]));
567 : cx[s] += (aux * cos(2.0 * M_PI * Position[0]));
568 : cy[s] += (aux * cos(2.0 * M_PI * Position[1]));
569 : }
570 : }
571 : }
572 : }
573 : }
574 : }
575 :
576 184 : for (unsigned s = 0; s < NS[0]; s++)
577 : {
578 180 : if (Fs[s] != 0.0)
579 : {
580 49 : ws[s] = tanh(Fs[s]);
581 49 : W += ws[s];
582 49 : sx[s] = sx[s] / Fs[s];
583 49 : sy[s] = sy[s] / Fs[s];
584 49 : cx[s] = cx[s] / Fs[s];
585 49 : cy[s] = cy[s] / Fs[s];
586 49 : Xsc += sx[s] * ws[s];
587 49 : Ysc += sy[s] * ws[s];
588 49 : Xcc += cx[s] * ws[s];
589 49 : Ycc += cy[s] * ws[s];
590 : }
591 : }
592 :
593 4 : Xsc = Xsc / W;
594 4 : Ysc = Ysc / W;
595 4 : Xcc = Xcc / W;
596 4 : Ycc = Ycc / W;
597 :
598 : // Eq. 12 Hub & Awasthi JCTC 2017.
599 : double Xcyl, Ycyl;
600 :
601 : Xcyl = Xcyl_Mem;
602 : Ycyl = Ycyl_Mem;
603 :
604 : // Eq. 25, 26 and 27 Hub & Awasthi JCTC 2017.
605 : double d_sx_dx, d_sx_dz, d_sy_dy, d_sy_dz, d_cx_dx, d_cx_dz, d_cy_dy, d_cy_dz;
606 :
607 : // Eq. 29 Hub & Awasthi JCTC 2017.
608 : double d_ws_dz;
609 :
610 : // Eq. 31, 32 and 33 Hub & Awasthi JCTC 2017
611 : double d_Xsc_dx, d_Xsc_dz, d_Xcc_dx, d_Xcc_dz, d_Ysc_dy, d_Ysc_dz, d_Ycc_dy, d_Ycc_dz;
612 :
613 : // Center of the cylinder. X and Y are calculated (or defined), Z is the Z component of the geometric center of the membranes.
614 4 : Vector xyzCyl = pbcDistance(Vector(0.0, 0.0, 0.0), Vector(Xcyl, Ycyl, ZMems));
615 :
616 : // Distances from the oxygens to center of the cylinder.
617 4 : std::vector<Vector> CylDistances(chainBeads);
618 :
619 : // Modulo of the XY distances from the oxygens to the center of the cylinder.
620 : double ri;
621 :
622 : // Eq. 8 Hub & Awasthi JCTC 2017.
623 : double fradial;
624 :
625 : // Eq. 15 Hub & Awasthi JCTC 2017.
626 4 : std::vector<double> d_fradial_dx(chainBeads), d_fradial_dy(chainBeads);
627 :
628 : // Eq. 35, 36, 37 and 38 Hub & Awasthi JCTC 2017.
629 4 : std::vector<double> d_Xcyl_dx(chainBeads), d_Xcyl_dz(chainBeads), d_Ycyl_dy(chainBeads), d_Ycyl_dz(chainBeads);
630 :
631 : // To avoid rare instabilities auxX and auxY are truncated at a configurable value (default 500).
632 4 : double auxX = (1 / (pow(Xsc, 2) + pow(Xcc, 2))), auxY = (1 / (pow(Ysc, 2) + pow(Ycc, 2)));
633 :
634 4 : if (auxX > ONEOVERS2C2CUTOFF[0])
635 : {
636 0 : auxX = Lx * ONEOVERS2C2CUTOFF[0] / (2 * M_PI);
637 : }
638 : else
639 : {
640 4 : auxX = Lx * auxX / (2 * M_PI);
641 : }
642 :
643 4 : if (auxY > ONEOVERS2C2CUTOFF[0])
644 : {
645 0 : auxY = Ly * ONEOVERS2C2CUTOFF[0] / (2 * M_PI);
646 : }
647 : else
648 : {
649 4 : auxY = Ly * auxY / (2 * M_PI);
650 : }
651 :
652 : //Number of oxygens within the slice s of the membrane-spanning cylinder.
653 4 : std::vector<double> Nsp(NS[0]), psi(NS[0]), d_psi(NS[0]);
654 :
655 : // Eq. 3 Hub & Awasthi JCTC 2017.
656 4 : double b = (ZETA[0] / (1.0 - ZETA[0])), c = ((1.0 - ZETA[0]) * exp(b));
657 :
658 : // Eq. 19 Hub & Awasthi JCTC 2017.
659 4 : std::vector<double> fradial_d_faxial_dz(chainBeads * NS[0]);
660 :
661 : // Eq. 20 Hub & Awasthi JCTC 2017.
662 4 : std::vector<double> Axs(NS[0]), Ays(NS[0]);
663 :
664 : #ifdef _OPENMP
665 : #if _OPENMP >= 201307
666 4 : #pragma omp parallel for private(d_Xsc_dx,d_Xcc_dx,d_Ysc_dy,d_Ycc_dy,d_Xsc_dz,d_Xcc_dz,d_Ysc_dz,d_Ycc_dz,d_sx_dx,d_sy_dy,d_cx_dx,d_cy_dy,d_sx_dz,d_sy_dz,d_cx_dz,d_cy_dz,d_ws_dz,ri,x,fradial) reduction(vec_double_plus: Nsp, Axs, Ays)
667 : #endif
668 : #endif
669 : for (unsigned i = 0; i < chainBeads; i++)
670 : {
671 : CylDistances[i] = pbcDistance(xyzCyl, pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + noChainBeads)));
672 : if (analyzeThisParticle[i])
673 : {
674 : Position = getPbc().realToScaled(pbcDistance(Vector(0.0, 0.0, 0.0), getPosition(i + noChainBeads)));
675 : d_Xsc_dx = 0.0;
676 : d_Xcc_dx = 0.0;
677 : d_Ysc_dy = 0.0;
678 : d_Ycc_dy = 0.0;
679 : d_Xsc_dz = 0.0;
680 : d_Xcc_dz = 0.0;
681 : d_Ysc_dz = 0.0;
682 : d_Ycc_dz = 0.0;
683 : for (unsigned s = s1[i]; s <= s2[i]; s++)
684 : {
685 : if (Fs[s] != 0.0)
686 : {
687 : d_sx_dx = faxial[i + chainBeads * s] * 2.0 * M_PI * cos(2.0 * M_PI * Position[0]) / (Lx * Fs[s]);
688 : d_sy_dy = faxial[i + chainBeads * s] * 2.0 * M_PI * cos(2.0 * M_PI * Position[1]) / (Ly * Fs[s]);
689 : d_cx_dx = -faxial[i + chainBeads * s] * 2.0 * M_PI * sin(2.0 * M_PI * Position[0]) / (Lx * Fs[s]);
690 : d_cy_dy = -faxial[i + chainBeads * s] * 2.0 * M_PI * sin(2.0 * M_PI * Position[1]) / (Ly * Fs[s]);
691 : d_Xsc_dx += ws[s] * d_sx_dx / W;
692 : d_Xcc_dx += ws[s] * d_cx_dx / W;
693 : d_Ysc_dy += ws[s] * d_sy_dy / W;
694 : d_Ycc_dy += ws[s] * d_cy_dy / W;
695 :
696 : d_sx_dz = d_faxial_dz[i + chainBeads * s] * (sin(2.0 * M_PI * Position[0]) - sx[s]) / Fs[s];
697 : d_sy_dz = d_faxial_dz[i + chainBeads * s] * (sin(2.0 * M_PI * Position[1]) - sy[s]) / Fs[s];
698 : d_cx_dz = d_faxial_dz[i + chainBeads * s] * (cos(2.0 * M_PI * Position[0]) - cx[s]) / Fs[s];
699 : d_cy_dz = d_faxial_dz[i + chainBeads * s] * (cos(2.0 * M_PI * Position[1]) - cy[s]) / Fs[s];
700 : d_ws_dz = (1 - pow(ws[s], 2)) * d_faxial_dz[i + chainBeads * s];
701 : d_Xsc_dz += (ws[s] * d_sx_dz + d_ws_dz * (sx[s] - Xsc)) / W;
702 : d_Xcc_dz += (ws[s] * d_cx_dz + d_ws_dz * (cx[s] - Xcc)) / W;
703 : d_Ysc_dz += (ws[s] * d_sy_dz + d_ws_dz * (sy[s] - Ysc)) / W;
704 : d_Ycc_dz += (ws[s] * d_cy_dz + d_ws_dz * (cy[s] - Ycc)) / W;
705 : }
706 : }
707 : d_Xcyl_dx[i] = auxX * (-Xsc * d_Xcc_dx + Xcc * d_Xsc_dx);
708 : d_Xcyl_dz[i] = auxX * (-Xsc * d_Xcc_dz + Xcc * d_Xsc_dz);
709 : d_Ycyl_dy[i] = auxY * (-Ysc * d_Ycc_dy + Ycc * d_Ysc_dy);
710 : d_Ycyl_dz[i] = auxY * (-Ysc * d_Ycc_dz + Ycc * d_Ysc_dz);
711 :
712 : ri = sqrt(pow(CylDistances[i][0], 2) + pow(CylDistances[i][1], 2));
713 : x = ri / RCYL[0];
714 : if (!((x <= -1.0 - HCH[0]) || (x >= 1.0 + HCH[0])))
715 : {
716 : if (((-1.0 + HCH[0]) <= x) && (x <= (1.0 - HCH[0])))
717 : {
718 : fradial = 1.0;
719 : }
720 : else if (((1.0 - HCH[0]) < x) && (x < (1.0 + HCH[0])))
721 : {
722 : fradial = 0.5 - ((3.0 * x - 3.0) / (4.0 * HCH[0])) + (pow((x - 1.0), 3) / (4.0 * pow(HCH[0], 3)));
723 : d_fradial_dx[i] = ((-3.0 / (4.0 * HCH[0])) + ((3.0 * pow((x - 1), 2)) / (4.0 * pow(HCH[0], 3)))) * CylDistances[i][0] / (RCYL[0] * ri);
724 : d_fradial_dy[i] = ((-3.0 / (4.0 * HCH[0])) + ((3.0 * pow((x - 1), 2)) / (4.0 * pow(HCH[0], 3)))) * CylDistances[i][1] / (RCYL[0] * ri);
725 : }
726 : else if (((-1.0 - HCH[0]) < x) && (x < (-1.0 + HCH[0])))
727 : {
728 : fradial = 0.5 + ((3.0 * x + 3.0) / (4.0 * HCH[0])) - (pow((x + 1.0), 3) / (4.0 * pow(HCH[0], 3)));
729 : d_fradial_dx[i] = ((3.0 / (4.0 * HCH[0])) - ((3.0 * pow((x + 1), 2)) / (4.0 * pow(HCH[0], 3)))) * CylDistances[i][0] / (RCYL[0] * ri);
730 : d_fradial_dy[i] = ((3.0 / (4.0 * HCH[0])) - ((3.0 * pow((x + 1), 2)) / (4.0 * pow(HCH[0], 3)))) * CylDistances[i][1] / (RCYL[0] * ri);
731 : }
732 :
733 : for (unsigned s = s1[i]; s <= s2[i]; s++)
734 : {
735 : Nsp[s] += fradial * faxial[i + chainBeads * s];
736 : Axs[s] += faxial[i + chainBeads * s] * d_fradial_dx[i];
737 : Ays[s] += faxial[i + chainBeads * s] * d_fradial_dy[i];
738 : fradial_d_faxial_dz[i + chainBeads * s] = fradial * d_faxial_dz[i + chainBeads * s];
739 : }
740 : }
741 : }
742 : }
743 :
744 184 : for (unsigned s = 0; s < NS[0]; s++)
745 : {
746 180 : if (Nsp[s] <= 1.0)
747 : {
748 149 : psi[s] = ZETA[0] * Nsp[s];
749 149 : d_psi[s] = ZETA[0];
750 149 : Xi_n += psi[s];
751 : }
752 : else
753 : {
754 31 : psi[s] = 1.0 - c * exp(-b * Nsp[s]);
755 31 : d_psi[s] = b * c * exp(-b * Nsp[s]);
756 31 : Xi_n += psi[s];
757 : }
758 : }
759 :
760 4 : Xi_n = Xi_n / NS[0];
761 :
762 : // Eq. 18 Hub & Awasthi JCTC 2017.
763 4 : std::vector<double> faxial_d_fradial_dx(chainBeads * NS[0]), faxial_d_fradial_dy(chainBeads * NS[0]), faxial_d_fradial_dz(chainBeads * NS[0]);
764 :
765 : // Eq. 13 Hub & Awasthi JCTC 2017 modified to considere the Heaviside_Chain step function (this only affect during the transition).
766 4 : std::vector<Vector> derivatives_Chain(chainBeads);
767 :
768 : #ifdef _OPENMP
769 : #if _OPENMP >= 201307
770 4 : #pragma omp parallel for private(aux)
771 : #endif
772 : #endif
773 : for (unsigned i = 0; i < chainBeads; i++)
774 : {
775 : if (analyzeThisParticle[i])
776 : {
777 : for (unsigned s = s1[i]; s <= s2[i]; s++)
778 : {
779 : if (faxial[i + chainBeads * s])
780 : {
781 : faxial_d_fradial_dx[i + chainBeads * s] = faxial[i + chainBeads * s] * d_fradial_dx[i] - d_Xcyl_dx[i] * Axs[s];
782 : faxial_d_fradial_dy[i + chainBeads * s] = faxial[i + chainBeads * s] * d_fradial_dy[i] - d_Ycyl_dy[i] * Ays[s];
783 : faxial_d_fradial_dz[i + chainBeads * s] = -d_Xcyl_dz[i] * Axs[s] - d_Ycyl_dz[i] * Ays[s];
784 : }
785 : }
786 :
787 : for (unsigned s = s1[i]; s <= s2[i]; s++)
788 : {
789 : aux = d_psi[s] / NS[0];
790 : derivatives_Chain[i][0] += aux * faxial_d_fradial_dx[i + chainBeads * s];
791 : derivatives_Chain[i][1] += aux * faxial_d_fradial_dy[i + chainBeads * s];
792 : derivatives_Chain[i][2] += aux * (faxial_d_fradial_dz[i + chainBeads * s] + fradial_d_faxial_dz[i + chainBeads * s]);
793 : }
794 : }
795 : }
796 :
797 4 : Tensor virial;
798 134604 : for (unsigned i = 0; i < chainBeads; i++)
799 : {
800 134600 : setAtomsDerivatives((i + noChainBeads), derivatives_Chain[i]);
801 134600 : virial -= Tensor(CylDistances[i], derivatives_Chain[i]);
802 : }
803 :
804 4 : setValue(Xi_n);
805 4 : setBoxDerivatives(virial);
806 4 : }
807 : }
808 : }
|
