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