2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2008, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
49 #include "gromacs/domdec/dlbtiming.h"
50 #include "gromacs/domdec/domdec.h"
51 #include "gromacs/domdec/domdec_struct.h"
52 #include "gromacs/gmxlib/chargegroup.h"
53 #include "gromacs/gmxlib/network.h"
54 #include "gromacs/math/functions.h"
55 #include "gromacs/math/units.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/math/vecdump.h"
58 #include "gromacs/mdlib/constr.h"
59 #include "gromacs/mdlib/force.h"
60 #include "gromacs/mdlib/mdrun.h"
61 #include "gromacs/mdlib/sim_util.h"
62 #include "gromacs/mdlib/vsite.h"
63 #include "gromacs/mdtypes/commrec.h"
64 #include "gromacs/mdtypes/inputrec.h"
65 #include "gromacs/mdtypes/md_enums.h"
66 #include "gromacs/mdtypes/state.h"
67 #include "gromacs/pbcutil/mshift.h"
68 #include "gromacs/pbcutil/pbc.h"
69 #include "gromacs/topology/ifunc.h"
70 #include "gromacs/topology/mtop_lookup.h"
71 #include "gromacs/topology/mtop_util.h"
72 #include "gromacs/utility/arrayref.h"
73 #include "gromacs/utility/arraysize.h"
74 #include "gromacs/utility/cstringutil.h"
75 #include "gromacs/utility/fatalerror.h"
76 #include "gromacs/utility/smalloc.h"
80 int shell
; /* The shell id */
81 int nucl1
, nucl2
, nucl3
; /* The nuclei connected to the shell */
82 /* gmx_bool bInterCG; */ /* Coupled to nuclei outside cg? */
83 real k
; /* force constant */
84 real k_1
; /* 1 over force constant */
90 struct gmx_shellfc_t
{
91 /* Shell counts, indices, parameters and working data */
92 int nshell_gl
; /* The number of shells in the system */
93 t_shell
*shell_gl
; /* All the shells (for DD only) */
94 int *shell_index_gl
; /* Global shell index (for DD only) */
95 gmx_bool bInterCG
; /* Are there inter charge-group shells? */
96 int nshell
; /* The number of local shells */
97 t_shell
*shell
; /* The local shells */
98 int shell_nalloc
; /* The allocation size of shell */
99 gmx_bool bPredict
; /* Predict shell positions */
100 gmx_bool bRequireInit
; /* Require initialization of shell positions */
101 int nflexcon
; /* The number of flexible constraints */
103 /* Temporary arrays, should be fixed size 2 when fully converted to C++ */
104 PaddedRVecVector
*x
; /* Array for iterative minimization */
105 PaddedRVecVector
*f
; /* Array for iterative minimization */
107 /* Flexible constraint working data */
108 rvec
*acc_dir
; /* Acceleration direction for flexcon */
109 rvec
*x_old
; /* Old coordinates for flexcon */
110 int flex_nalloc
; /* The allocation size of acc_dir and x_old */
111 rvec
*adir_xnold
; /* Work space for init_adir */
112 rvec
*adir_xnew
; /* Work space for init_adir */
113 int adir_nalloc
; /* Work space for init_adir */
114 std::int64_t numForceEvaluations
; /* Total number of force evaluations */
115 int numConvergedIterations
; /* Total number of iterations that converged */
119 static void pr_shell(FILE *fplog
, int ns
, t_shell s
[])
123 fprintf(fplog
, "SHELL DATA\n");
124 fprintf(fplog
, "%5s %8s %5s %5s %5s\n",
125 "Shell", "Force k", "Nucl1", "Nucl2", "Nucl3");
126 for (i
= 0; (i
< ns
); i
++)
128 fprintf(fplog
, "%5d %8.3f %5d", s
[i
].shell
, 1.0/s
[i
].k_1
, s
[i
].nucl1
);
131 fprintf(fplog
, " %5d\n", s
[i
].nucl2
);
133 else if (s
[i
].nnucl
== 3)
135 fprintf(fplog
, " %5d %5d\n", s
[i
].nucl2
, s
[i
].nucl3
);
139 fprintf(fplog
, "\n");
144 /* TODO The remain call of this function passes non-NULL mass and NULL
145 * mtop, so this routine can be simplified.
147 * The other code path supported doing prediction before the MD loop
148 * started, but even when called, the prediction was always
149 * over-written by a subsequent call in the MD loop, so has been
151 static void predict_shells(FILE *fplog
, rvec x
[], rvec v
[], real dt
,
153 real mass
[], gmx_mtop_t
*mtop
, gmx_bool bInit
)
155 int i
, m
, s1
, n1
, n2
, n3
;
156 real dt_1
, fudge
, tm
, m1
, m2
, m3
;
159 /* We introduce a fudge factor for performance reasons: with this choice
160 * the initial force on the shells is about a factor of two lower than
169 fprintf(fplog
, "RELAX: Using prediction for initial shell placement\n");
181 for (i
= 0; (i
< ns
); i
++)
192 for (m
= 0; (m
< DIM
); m
++)
194 x
[s1
][m
] += ptr
[n1
][m
]*dt_1
;
207 /* Not the correct masses with FE, but it is just a prediction... */
208 m1
= mtopGetAtomMass(mtop
, n1
, &molb
);
209 m2
= mtopGetAtomMass(mtop
, n2
, &molb
);
212 for (m
= 0; (m
< DIM
); m
++)
214 x
[s1
][m
] += (m1
*ptr
[n1
][m
]+m2
*ptr
[n2
][m
])*tm
;
229 /* Not the correct masses with FE, but it is just a prediction... */
230 m1
= mtopGetAtomMass(mtop
, n1
, &molb
);
231 m2
= mtopGetAtomMass(mtop
, n2
, &molb
);
232 m3
= mtopGetAtomMass(mtop
, n3
, &molb
);
234 tm
= dt_1
/(m1
+m2
+m3
);
235 for (m
= 0; (m
< DIM
); m
++)
237 x
[s1
][m
] += (m1
*ptr
[n1
][m
]+m2
*ptr
[n2
][m
]+m3
*ptr
[n3
][m
])*tm
;
241 gmx_fatal(FARGS
, "Shell %d has %d nuclei!", i
, s
[i
].nnucl
);
246 /*! \brief Count the different particle types in a system
248 * Routine prints a warning to stderr in case an unknown particle type
250 * \param[in] fplog Print what we have found if not NULL
251 * \param[in] mtop Molecular topology.
252 * \returns Array holding the number of particles of a type
254 static std::array
<int, eptNR
> countPtypes(FILE *fplog
,
255 const gmx_mtop_t
*mtop
)
257 std::array
<int, eptNR
> nptype
= { { 0 } };
258 /* Count number of shells, and find their indices */
259 for (int i
= 0; (i
< eptNR
); i
++)
264 gmx_mtop_atomloop_block_t aloopb
= gmx_mtop_atomloop_block_init(mtop
);
267 while (gmx_mtop_atomloop_block_next(aloopb
, &atom
, &nmol
))
274 nptype
[atom
->ptype
] += nmol
;
277 fprintf(stderr
, "Warning unsupported particle type %d in countPtypes",
278 static_cast<int>(atom
->ptype
));
283 /* Print the number of each particle type */
285 for (const auto &i
: nptype
)
289 fprintf(fplog
, "There are: %d %ss\n", i
, ptype_str
[n
]);
297 gmx_shellfc_t
*init_shell_flexcon(FILE *fplog
,
298 const gmx_mtop_t
*mtop
, int nflexcon
,
300 bool usingDomainDecomposition
)
304 int *shell_index
= nullptr, *at2cg
;
308 int i
, j
, type
, a_offset
, cg
, mol
, ftype
, nra
;
310 int aS
, aN
= 0; /* Shell and nucleus */
311 int bondtypes
[] = { F_BONDS
, F_HARMONIC
, F_CUBICBONDS
, F_POLARIZATION
, F_ANHARM_POL
, F_WATER_POL
};
312 #define NBT asize(bondtypes)
314 gmx_mtop_atomloop_all_t aloop
;
315 const gmx_ffparams_t
*ffparams
;
317 std::array
<int, eptNR
> n
= countPtypes(fplog
, mtop
);
318 nshell
= n
[eptShell
];
320 if (nshell
== 0 && nflexcon
== 0)
322 /* We're not doing shells or flexible constraints */
327 shfc
->x
= new PaddedRVecVector
[2] {};
328 shfc
->f
= new PaddedRVecVector
[2] {};
329 shfc
->nflexcon
= nflexcon
;
333 /* Only flexible constraints, no shells.
334 * Note that make_local_shells() does not need to be called.
337 shfc
->bPredict
= FALSE
;
342 if (nstcalcenergy
!= 1)
344 gmx_fatal(FARGS
, "You have nstcalcenergy set to a value (%d) that is different from 1.\nThis is not supported in combination with shell particles.\nPlease make a new tpr file.", nstcalcenergy
);
346 if (usingDomainDecomposition
)
348 gmx_fatal(FARGS
, "Shell particles are not implemented with domain decomposition, use a single rank");
351 /* We have shells: fill the shell data structure */
353 /* Global system sized array, this should be avoided */
354 snew(shell_index
, mtop
->natoms
);
356 aloop
= gmx_mtop_atomloop_all_init(mtop
);
358 while (gmx_mtop_atomloop_all_next(aloop
, &i
, &atom
))
360 if (atom
->ptype
== eptShell
)
362 shell_index
[i
] = nshell
++;
368 /* Initiate the shell structures */
369 for (i
= 0; (i
< nshell
); i
++)
376 /* shell[i].bInterCG=FALSE; */
381 ffparams
= &mtop
->ffparams
;
383 /* Now fill the structures */
384 shfc
->bInterCG
= FALSE
;
387 for (size_t mb
= 0; mb
< mtop
->molblock
.size(); mb
++)
389 const gmx_molblock_t
*molb
= &mtop
->molblock
[mb
];
390 const gmx_moltype_t
*molt
= &mtop
->moltype
[molb
->type
];
391 const t_block
*cgs
= &molt
->cgs
;
393 snew(at2cg
, molt
->atoms
.nr
);
394 for (cg
= 0; cg
< cgs
->nr
; cg
++)
396 for (i
= cgs
->index
[cg
]; i
< cgs
->index
[cg
+1]; i
++)
402 atom
= molt
->atoms
.atom
;
403 for (mol
= 0; mol
< molb
->nmol
; mol
++)
405 for (j
= 0; (j
< NBT
); j
++)
407 ia
= molt
->ilist
[bondtypes
[j
]].iatoms
;
408 for (i
= 0; (i
< molt
->ilist
[bondtypes
[j
]].nr
); )
411 ftype
= ffparams
->functype
[type
];
412 nra
= interaction_function
[ftype
].nratoms
;
414 /* Check whether we have a bond with a shell */
417 switch (bondtypes
[j
])
424 if (atom
[ia
[1]].ptype
== eptShell
)
429 else if (atom
[ia
[2]].ptype
== eptShell
)
436 aN
= ia
[4]; /* Dummy */
437 aS
= ia
[5]; /* Shell */
440 gmx_fatal(FARGS
, "Death Horror: %s, %d", __FILE__
, __LINE__
);
447 /* Check whether one of the particles is a shell... */
448 nsi
= shell_index
[a_offset
+aS
];
449 if ((nsi
< 0) || (nsi
>= nshell
))
451 gmx_fatal(FARGS
, "nsi is %d should be within 0 - %d. aS = %d",
454 if (shell
[nsi
].shell
== -1)
456 shell
[nsi
].shell
= a_offset
+ aS
;
459 else if (shell
[nsi
].shell
!= a_offset
+aS
)
461 gmx_fatal(FARGS
, "Weird stuff in %s, %d", __FILE__
, __LINE__
);
464 if (shell
[nsi
].nucl1
== -1)
466 shell
[nsi
].nucl1
= a_offset
+ aN
;
468 else if (shell
[nsi
].nucl2
== -1)
470 shell
[nsi
].nucl2
= a_offset
+ aN
;
472 else if (shell
[nsi
].nucl3
== -1)
474 shell
[nsi
].nucl3
= a_offset
+ aN
;
480 pr_shell(fplog
, ns
, shell
);
482 gmx_fatal(FARGS
, "Can not handle more than three bonds per shell\n");
484 if (at2cg
[aS
] != at2cg
[aN
])
486 /* shell[nsi].bInterCG = TRUE; */
487 shfc
->bInterCG
= TRUE
;
490 switch (bondtypes
[j
])
494 shell
[nsi
].k
+= ffparams
->iparams
[type
].harmonic
.krA
;
497 shell
[nsi
].k
+= ffparams
->iparams
[type
].cubic
.kb
;
501 if (!gmx_within_tol(qS
, atom
[aS
].qB
, GMX_REAL_EPS
*10))
503 gmx_fatal(FARGS
, "polarize can not be used with qA(%e) != qB(%e) for atom %d of molecule block %d", qS
, atom
[aS
].qB
, aS
+1, mb
+1);
505 shell
[nsi
].k
+= gmx::square(qS
)*ONE_4PI_EPS0
/
506 ffparams
->iparams
[type
].polarize
.alpha
;
509 if (!gmx_within_tol(qS
, atom
[aS
].qB
, GMX_REAL_EPS
*10))
511 gmx_fatal(FARGS
, "water_pol can not be used with qA(%e) != qB(%e) for atom %d of molecule block %d", qS
, atom
[aS
].qB
, aS
+1, mb
+1);
513 alpha
= (ffparams
->iparams
[type
].wpol
.al_x
+
514 ffparams
->iparams
[type
].wpol
.al_y
+
515 ffparams
->iparams
[type
].wpol
.al_z
)/3.0;
516 shell
[nsi
].k
+= gmx::square(qS
)*ONE_4PI_EPS0
/alpha
;
519 gmx_fatal(FARGS
, "Death Horror: %s, %d", __FILE__
, __LINE__
);
527 a_offset
+= molt
->atoms
.nr
;
529 /* Done with this molecule type */
533 /* Verify whether it's all correct */
536 gmx_fatal(FARGS
, "Something weird with shells. They may not be bonded to something");
539 for (i
= 0; (i
< ns
); i
++)
541 shell
[i
].k_1
= 1.0/shell
[i
].k
;
546 pr_shell(debug
, ns
, shell
);
550 shfc
->nshell_gl
= ns
;
551 shfc
->shell_gl
= shell
;
552 shfc
->shell_index_gl
= shell_index
;
554 shfc
->bPredict
= (getenv("GMX_NOPREDICT") == nullptr);
555 shfc
->bRequireInit
= FALSE
;
560 fprintf(fplog
, "\nWill never predict shell positions\n");
565 shfc
->bRequireInit
= (getenv("GMX_REQUIRE_SHELL_INIT") != nullptr);
566 if (shfc
->bRequireInit
&& fplog
)
568 fprintf(fplog
, "\nWill always initiate shell positions\n");
578 fprintf(fplog
, "\nNOTE: there all shells that are connected to particles outside thier own charge group, will not predict shells positions during the run\n\n");
580 /* Prediction improves performance, so we should implement either:
581 * 1. communication for the atoms needed for prediction
582 * 2. prediction using the velocities of shells; currently the
583 * shell velocities are zeroed, it's a bit tricky to keep
584 * track of the shell displacements and thus the velocity.
586 shfc
->bPredict
= FALSE
;
593 void make_local_shells(const t_commrec
*cr
, t_mdatoms
*md
,
597 int a0
, a1
, *ind
, nshell
, i
;
598 gmx_domdec_t
*dd
= nullptr;
600 if (DOMAINDECOMP(cr
))
608 /* Single node: we need all shells, just copy the pointer */
609 shfc
->nshell
= shfc
->nshell_gl
;
610 shfc
->shell
= shfc
->shell_gl
;
615 ind
= shfc
->shell_index_gl
;
619 for (i
= a0
; i
< a1
; i
++)
621 if (md
->ptype
[i
] == eptShell
)
623 if (nshell
+1 > shfc
->shell_nalloc
)
625 shfc
->shell_nalloc
= over_alloc_dd(nshell
+1);
626 srenew(shell
, shfc
->shell_nalloc
);
630 shell
[nshell
] = shfc
->shell_gl
[ind
[dd
->gatindex
[i
]]];
634 shell
[nshell
] = shfc
->shell_gl
[ind
[i
]];
637 /* With inter-cg shells we can no do shell prediction,
638 * so we do not need the nuclei numbers.
642 shell
[nshell
].nucl1
= i
+ shell
[nshell
].nucl1
- shell
[nshell
].shell
;
643 if (shell
[nshell
].nnucl
> 1)
645 shell
[nshell
].nucl2
= i
+ shell
[nshell
].nucl2
- shell
[nshell
].shell
;
647 if (shell
[nshell
].nnucl
> 2)
649 shell
[nshell
].nucl3
= i
+ shell
[nshell
].nucl3
- shell
[nshell
].shell
;
652 shell
[nshell
].shell
= i
;
657 shfc
->nshell
= nshell
;
661 static void do_1pos(rvec xnew
, const rvec xold
, const rvec f
, real step
)
679 static void do_1pos3(rvec xnew
, const rvec xold
, const rvec f
, const rvec step
)
697 static void directional_sd(gmx::ArrayRef
<const gmx::RVec
> xold
,
698 gmx::ArrayRef
<gmx::RVec
> xnew
,
699 const rvec acc_dir
[], int homenr
, real step
)
701 const rvec
*xo
= as_rvec_array(xold
.data());
702 rvec
*xn
= as_rvec_array(xnew
.data());
704 for (int i
= 0; i
< homenr
; i
++)
706 do_1pos(xn
[i
], xo
[i
], acc_dir
[i
], step
);
710 static void shell_pos_sd(gmx::ArrayRef
<const gmx::RVec
> xcur
,
711 gmx::ArrayRef
<gmx::RVec
> xnew
,
712 gmx::ArrayRef
<const gmx::RVec
> f
,
713 int ns
, t_shell s
[], int count
)
715 const real step_scale_min
= 0.8,
716 step_scale_increment
= 0.2,
717 step_scale_max
= 1.2,
718 step_scale_multiple
= (step_scale_max
- step_scale_min
) / step_scale_increment
;
723 real step_min
, step_max
;
728 for (i
= 0; (i
< ns
); i
++)
733 for (d
= 0; d
< DIM
; d
++)
735 s
[i
].step
[d
] = s
[i
].k_1
;
737 step_min
= std::min(step_min
, s
[i
].step
[d
]);
738 step_max
= std::max(step_max
, s
[i
].step
[d
]);
744 for (d
= 0; d
< DIM
; d
++)
746 dx
= xcur
[shell
][d
] - s
[i
].xold
[d
];
747 df
= f
[shell
][d
] - s
[i
].fold
[d
];
748 /* -dx/df gets used to generate an interpolated value, but would
749 * cause a NaN if df were binary-equal to zero. Values close to
750 * zero won't cause problems (because of the min() and max()), so
751 * just testing for binary inequality is OK. */
755 /* Scale the step size by a factor interpolated from
756 * step_scale_min to step_scale_max, as k_est goes from 0 to
757 * step_scale_multiple * s[i].step[d] */
759 step_scale_min
* s
[i
].step
[d
] +
760 step_scale_increment
* std::min(step_scale_multiple
* s
[i
].step
[d
], std::max(k_est
, zero
));
765 if (gmx_numzero(dx
)) /* 0 == dx */
767 /* Likely this will never happen, but if it does just
768 * don't scale the step. */
772 s
[i
].step
[d
] *= step_scale_max
;
776 step_min
= std::min(step_min
, s
[i
].step
[d
]);
777 step_max
= std::max(step_max
, s
[i
].step
[d
]);
781 copy_rvec(xcur
[shell
], s
[i
].xold
);
782 copy_rvec(f
[shell
], s
[i
].fold
);
784 do_1pos3(xnew
[shell
], xcur
[shell
], f
[shell
], s
[i
].step
);
788 fprintf(debug
, "shell[%d] = %d\n", i
, shell
);
789 pr_rvec(debug
, 0, "fshell", f
[shell
], DIM
, TRUE
);
790 pr_rvec(debug
, 0, "xold", xcur
[shell
], DIM
, TRUE
);
791 pr_rvec(debug
, 0, "step", s
[i
].step
, DIM
, TRUE
);
792 pr_rvec(debug
, 0, "xnew", xnew
[shell
], DIM
, TRUE
);
796 printf("step %.3e %.3e\n", step_min
, step_max
);
800 static void decrease_step_size(int nshell
, t_shell s
[])
804 for (i
= 0; i
< nshell
; i
++)
806 svmul(0.8, s
[i
].step
, s
[i
].step
);
810 static void print_epot(FILE *fp
, gmx_int64_t mdstep
, int count
, real epot
, real df
,
811 int ndir
, real sf_dir
)
815 fprintf(fp
, "MDStep=%5s/%2d EPot: %12.8e, rmsF: %6.2e",
816 gmx_step_str(mdstep
, buf
), count
, epot
, df
);
819 fprintf(fp
, ", dir. rmsF: %6.2e\n", std::sqrt(sf_dir
/ndir
));
828 static real
rms_force(const t_commrec
*cr
, gmx::ArrayRef
<const gmx::RVec
> force
, int ns
, t_shell s
[],
829 int ndir
, real
*sf_dir
, real
*Epot
)
832 const rvec
*f
= as_rvec_array(force
.data());
835 for (int i
= 0; i
< ns
; i
++)
837 int shell
= s
[i
].shell
;
838 buf
[0] += norm2(f
[shell
]);
847 gmx_sumd(4, buf
, cr
);
848 ntot
= (int)(buf
[1] + 0.5);
854 return (ntot
? std::sqrt(buf
[0]/ntot
) : 0);
857 static void check_pbc(FILE *fp
, gmx::ArrayRef
<gmx::RVec
> x
, int shell
)
862 for (m
= 0; (m
< DIM
); m
++)
864 if (fabs(x
[shell
][m
]-x
[now
][m
]) > 0.3)
866 pr_rvecs(fp
, 0, "SHELL-X", as_rvec_array(x
.data())+now
, 5);
872 static void dump_shells(FILE *fp
, gmx::ArrayRef
<gmx::RVec
> x
, gmx::ArrayRef
<gmx::RVec
> f
, real ftol
, int ns
, t_shell s
[])
877 ft2
= gmx::square(ftol
);
879 for (i
= 0; (i
< ns
); i
++)
882 ff2
= iprod(f
[shell
], f
[shell
]);
885 fprintf(fp
, "SHELL %5d, force %10.5f %10.5f %10.5f, |f| %10.5f\n",
886 shell
, f
[shell
][XX
], f
[shell
][YY
], f
[shell
][ZZ
], std::sqrt(ff2
));
888 check_pbc(fp
, x
, shell
);
892 static void init_adir(FILE *log
, gmx_shellfc_t
*shfc
,
893 gmx_constr_t constr
, t_idef
*idef
, t_inputrec
*ir
,
895 const gmx_multisim_t
*ms
,
897 gmx_int64_t step
, t_mdatoms
*md
, int end
,
898 rvec
*x_old
, rvec
*x_init
, rvec
*x
,
899 rvec
*f
, rvec
*acc_dir
,
900 gmx_bool bMolPBC
, matrix box
,
901 gmx::ArrayRef
<const real
> lambda
, real
*dvdlambda
,
907 unsigned short *ptype
;
909 if (DOMAINDECOMP(cr
))
917 if (n
> shfc
->adir_nalloc
)
919 shfc
->adir_nalloc
= over_alloc_dd(n
);
920 srenew(shfc
->adir_xnold
, shfc
->adir_nalloc
);
921 srenew(shfc
->adir_xnew
, shfc
->adir_nalloc
);
923 xnold
= shfc
->adir_xnold
;
924 xnew
= shfc
->adir_xnew
;
930 /* Does NOT work with freeze or acceleration groups (yet) */
931 for (n
= 0; n
< end
; n
++)
933 w_dt
= md
->invmass
[n
]*dt
;
935 for (d
= 0; d
< DIM
; d
++)
937 if ((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
))
939 xnold
[n
][d
] = x
[n
][d
] - (x_init
[n
][d
] - x_old
[n
][d
]);
940 xnew
[n
][d
] = 2*x
[n
][d
] - x_old
[n
][d
] + f
[n
][d
]*w_dt
*dt
;
944 xnold
[n
][d
] = x
[n
][d
];
945 xnew
[n
][d
] = x
[n
][d
];
949 constrain(log
, FALSE
, FALSE
, constr
, idef
, ir
, cr
, ms
, step
, 0, 1.0, md
,
950 x
, xnold
, nullptr, bMolPBC
, box
,
951 lambda
[efptBONDED
], &(dvdlambda
[efptBONDED
]),
952 nullptr, nullptr, nrnb
, econqCoord
);
953 constrain(log
, FALSE
, FALSE
, constr
, idef
, ir
, cr
, ms
, step
, 0, 1.0, md
,
954 x
, xnew
, nullptr, bMolPBC
, box
,
955 lambda
[efptBONDED
], &(dvdlambda
[efptBONDED
]),
956 nullptr, nullptr, nrnb
, econqCoord
);
958 for (n
= 0; n
< end
; n
++)
960 for (d
= 0; d
< DIM
; d
++)
963 -(2*x
[n
][d
]-xnold
[n
][d
]-xnew
[n
][d
])/gmx::square(dt
)
964 - f
[n
][d
]*md
->invmass
[n
];
966 clear_rvec(acc_dir
[n
]);
969 /* Project the acceleration on the old bond directions */
970 constrain(log
, FALSE
, FALSE
, constr
, idef
, ir
, cr
, ms
, step
, 0, 1.0, md
,
971 x_old
, xnew
, acc_dir
, bMolPBC
, box
,
972 lambda
[efptBONDED
], &(dvdlambda
[efptBONDED
]),
973 nullptr, nullptr, nrnb
, econqDeriv_FlexCon
);
976 void relax_shell_flexcon(FILE *fplog
, const t_commrec
*cr
,
977 const gmx_multisim_t
*ms
,
979 gmx_int64_t mdstep
, t_inputrec
*inputrec
,
980 gmx_bool bDoNS
, int force_flags
,
983 gmx_enerdata_t
*enerd
, t_fcdata
*fcd
,
984 t_state
*state
, PaddedRVecVector
*f
,
987 t_nrnb
*nrnb
, gmx_wallcycle_t wcycle
,
989 gmx_groups_t
*groups
,
992 double t
, rvec mu_tot
,
994 DdOpenBalanceRegionBeforeForceComputation ddOpenBalanceRegion
,
995 DdCloseBalanceRegionAfterForceComputation ddCloseBalanceRegion
)
1000 rvec
*acc_dir
= nullptr, *x_old
= nullptr;
1001 real Epot
[2], df
[2];
1005 gmx_bool bCont
, bInit
, bConverged
;
1006 int nat
, dd_ac0
, dd_ac1
= 0, i
;
1007 int homenr
= md
->homenr
, end
= homenr
, cg0
, cg1
;
1008 int nflexcon
, number_steps
, d
, Min
= 0, count
= 0;
1009 #define Try (1-Min) /* At start Try = 1 */
1011 bCont
= (mdstep
== inputrec
->init_step
) && inputrec
->bContinuation
;
1012 bInit
= (mdstep
== inputrec
->init_step
) || shfc
->bRequireInit
;
1013 ftol
= inputrec
->em_tol
;
1014 number_steps
= inputrec
->niter
;
1015 nshell
= shfc
->nshell
;
1016 shell
= shfc
->shell
;
1017 nflexcon
= shfc
->nflexcon
;
1021 if (DOMAINDECOMP(cr
))
1023 nat
= dd_natoms_vsite(cr
->dd
);
1026 dd_get_constraint_range(cr
->dd
, &dd_ac0
, &dd_ac1
);
1027 nat
= std::max(nat
, dd_ac1
);
1032 nat
= state
->natoms
;
1035 for (i
= 0; (i
< 2); i
++)
1037 shfc
->x
[i
].resize(gmx::paddedRVecVectorSize(nat
));
1038 shfc
->f
[i
].resize(gmx::paddedRVecVectorSize(nat
));
1041 /* Create views that we can swap */
1042 gmx::PaddedArrayRef
<gmx::RVec
> pos
[2];
1043 gmx::PaddedArrayRef
<gmx::RVec
> force
[2];
1044 for (i
= 0; (i
< 2); i
++)
1046 pos
[i
] = shfc
->x
[i
];
1047 force
[i
] = shfc
->f
[i
];
1050 if (bDoNS
&& inputrec
->ePBC
!= epbcNONE
&& !DOMAINDECOMP(cr
))
1052 /* This is the only time where the coordinates are used
1053 * before do_force is called, which normally puts all
1054 * charge groups in the box.
1056 if (inputrec
->cutoff_scheme
== ecutsVERLET
)
1058 auto xRef
= makeArrayRef(state
->x
);
1059 put_atoms_in_box_omp(fr
->ePBC
, state
->box
, xRef
.subArray(0, md
->homenr
));
1065 put_charge_groups_in_box(fplog
, cg0
, cg1
, fr
->ePBC
, state
->box
,
1066 &(top
->cgs
), as_rvec_array(state
->x
.data()), fr
->cg_cm
);
1071 mk_mshift(fplog
, graph
, fr
->ePBC
, state
->box
, as_rvec_array(state
->x
.data()));
1075 /* After this all coordinate arrays will contain whole charge groups */
1078 shift_self(graph
, state
->box
, as_rvec_array(state
->x
.data()));
1083 if (nat
> shfc
->flex_nalloc
)
1085 shfc
->flex_nalloc
= over_alloc_dd(nat
);
1086 srenew(shfc
->acc_dir
, shfc
->flex_nalloc
);
1087 srenew(shfc
->x_old
, shfc
->flex_nalloc
);
1089 acc_dir
= shfc
->acc_dir
;
1090 x_old
= shfc
->x_old
;
1091 for (i
= 0; i
< homenr
; i
++)
1093 for (d
= 0; d
< DIM
; d
++)
1096 state
->x
[i
][d
] - state
->v
[i
][d
]*inputrec
->delta_t
;
1101 /* Do a prediction of the shell positions, when appropriate.
1102 * Without velocities (EM, NM, BD) we only do initial prediction.
1104 if (shfc
->bPredict
&& !bCont
&& (EI_STATE_VELOCITY(inputrec
->eI
) || bInit
))
1106 predict_shells(fplog
, as_rvec_array(state
->x
.data()), as_rvec_array(state
->v
.data()), inputrec
->delta_t
, nshell
, shell
,
1107 md
->massT
, nullptr, bInit
);
1110 /* do_force expected the charge groups to be in the box */
1113 unshift_self(graph
, state
->box
, as_rvec_array(state
->x
.data()));
1116 /* Calculate the forces first time around */
1119 pr_rvecs(debug
, 0, "x b4 do_force", as_rvec_array(state
->x
.data()), homenr
);
1121 do_force(fplog
, cr
, ms
, inputrec
, mdstep
, nrnb
, wcycle
, top
, groups
,
1122 state
->box
, state
->x
, &state
->hist
,
1123 force
[Min
], force_vir
, md
, enerd
, fcd
,
1124 state
->lambda
, graph
,
1125 fr
, vsite
, mu_tot
, t
, nullptr,
1126 (bDoNS
? GMX_FORCE_NS
: 0) | force_flags
,
1127 ddOpenBalanceRegion
, ddCloseBalanceRegion
);
1132 init_adir(fplog
, shfc
,
1133 constr
, idef
, inputrec
, cr
, ms
, dd_ac1
, mdstep
, md
, end
,
1134 shfc
->x_old
, as_rvec_array(state
->x
.data()), as_rvec_array(state
->x
.data()), as_rvec_array(force
[Min
].data()),
1136 fr
->bMolPBC
, state
->box
, state
->lambda
, &dum
, nrnb
);
1138 for (i
= 0; i
< end
; i
++)
1140 sf_dir
+= md
->massT
[i
]*norm2(shfc
->acc_dir
[i
]);
1144 Epot
[Min
] = enerd
->term
[F_EPOT
];
1146 df
[Min
] = rms_force(cr
, shfc
->f
[Min
], nshell
, shell
, nflexcon
, &sf_dir
, &Epot
[Min
]);
1150 fprintf(debug
, "df = %g %g\n", df
[Min
], df
[Try
]);
1155 pr_rvecs(debug
, 0, "force0", as_rvec_array(force
[Min
].data()), md
->nr
);
1158 if (nshell
+nflexcon
> 0)
1160 /* Copy x to pos[Min] & pos[Try]: during minimization only the
1161 * shell positions are updated, therefore the other particles must
1164 pos
[Min
] = state
->x
;
1165 pos
[Try
] = state
->x
;
1168 if (bVerbose
&& MASTER(cr
))
1170 print_epot(stdout
, mdstep
, 0, Epot
[Min
], df
[Min
], nflexcon
, sf_dir
);
1175 fprintf(debug
, "%17s: %14.10e\n",
1176 interaction_function
[F_EKIN
].longname
, enerd
->term
[F_EKIN
]);
1177 fprintf(debug
, "%17s: %14.10e\n",
1178 interaction_function
[F_EPOT
].longname
, enerd
->term
[F_EPOT
]);
1179 fprintf(debug
, "%17s: %14.10e\n",
1180 interaction_function
[F_ETOT
].longname
, enerd
->term
[F_ETOT
]);
1181 fprintf(debug
, "SHELLSTEP %s\n", gmx_step_str(mdstep
, sbuf
));
1184 /* First check whether we should do shells, or whether the force is
1185 * low enough even without minimization.
1187 bConverged
= (df
[Min
] < ftol
);
1189 for (count
= 1; (!(bConverged
) && (count
< number_steps
)); count
++)
1193 construct_vsites(vsite
, as_rvec_array(pos
[Min
].data()),
1194 inputrec
->delta_t
, as_rvec_array(state
->v
.data()),
1195 idef
->iparams
, idef
->il
,
1196 fr
->ePBC
, fr
->bMolPBC
, cr
, state
->box
);
1201 init_adir(fplog
, shfc
,
1202 constr
, idef
, inputrec
, cr
, ms
, dd_ac1
, mdstep
, md
, end
,
1203 x_old
, as_rvec_array(state
->x
.data()), as_rvec_array(pos
[Min
].data()), as_rvec_array(force
[Min
].data()), acc_dir
,
1204 fr
->bMolPBC
, state
->box
, state
->lambda
, &dum
, nrnb
);
1206 directional_sd(pos
[Min
], pos
[Try
], acc_dir
, end
, fr
->fc_stepsize
);
1209 /* New positions, Steepest descent */
1210 shell_pos_sd(pos
[Min
], pos
[Try
], force
[Min
], nshell
, shell
, count
);
1212 /* do_force expected the charge groups to be in the box */
1215 unshift_self(graph
, state
->box
, as_rvec_array(pos
[Try
].data()));
1220 pr_rvecs(debug
, 0, "RELAX: pos[Min] ", as_rvec_array(pos
[Min
].data()), homenr
);
1221 pr_rvecs(debug
, 0, "RELAX: pos[Try] ", as_rvec_array(pos
[Try
].data()), homenr
);
1223 /* Try the new positions */
1224 do_force(fplog
, cr
, ms
, inputrec
, 1, nrnb
, wcycle
,
1225 top
, groups
, state
->box
, pos
[Try
], &state
->hist
,
1226 force
[Try
], force_vir
,
1227 md
, enerd
, fcd
, state
->lambda
, graph
,
1228 fr
, vsite
, mu_tot
, t
, nullptr,
1230 ddOpenBalanceRegion
, ddCloseBalanceRegion
);
1234 pr_rvecs(debug
, 0, "RELAX: force[Min]", as_rvec_array(force
[Min
].data()), homenr
);
1235 pr_rvecs(debug
, 0, "RELAX: force[Try]", as_rvec_array(force
[Try
].data()), homenr
);
1240 init_adir(fplog
, shfc
,
1241 constr
, idef
, inputrec
, cr
, ms
, dd_ac1
, mdstep
, md
, end
,
1242 x_old
, as_rvec_array(state
->x
.data()), as_rvec_array(pos
[Try
].data()), as_rvec_array(force
[Try
].data()), acc_dir
,
1243 fr
->bMolPBC
, state
->box
, state
->lambda
, &dum
, nrnb
);
1245 for (i
= 0; i
< end
; i
++)
1247 sf_dir
+= md
->massT
[i
]*norm2(acc_dir
[i
]);
1251 Epot
[Try
] = enerd
->term
[F_EPOT
];
1253 df
[Try
] = rms_force(cr
, force
[Try
], nshell
, shell
, nflexcon
, &sf_dir
, &Epot
[Try
]);
1257 fprintf(debug
, "df = %g %g\n", df
[Min
], df
[Try
]);
1264 pr_rvecs(debug
, 0, "F na do_force", as_rvec_array(force
[Try
].data()), homenr
);
1268 fprintf(debug
, "SHELL ITER %d\n", count
);
1269 dump_shells(debug
, pos
[Try
], force
[Try
], ftol
, nshell
, shell
);
1273 if (bVerbose
&& MASTER(cr
))
1275 print_epot(stdout
, mdstep
, count
, Epot
[Try
], df
[Try
], nflexcon
, sf_dir
);
1278 bConverged
= (df
[Try
] < ftol
);
1280 if ((df
[Try
] < df
[Min
]))
1284 fprintf(debug
, "Swapping Min and Try\n");
1288 /* Correct the velocities for the flexible constraints */
1289 invdt
= 1/inputrec
->delta_t
;
1290 for (i
= 0; i
< end
; i
++)
1292 for (d
= 0; d
< DIM
; d
++)
1294 state
->v
[i
][d
] += (pos
[Try
][i
][d
] - pos
[Min
][i
][d
])*invdt
;
1302 decrease_step_size(nshell
, shell
);
1305 shfc
->numForceEvaluations
+= count
;
1308 shfc
->numConvergedIterations
++;
1310 if (MASTER(cr
) && !(bConverged
))
1312 /* Note that the energies and virial are incorrect when not converged */
1316 "step %s: EM did not converge in %d iterations, RMS force %.3f\n",
1317 gmx_step_str(mdstep
, sbuf
), number_steps
, df
[Min
]);
1320 "step %s: EM did not converge in %d iterations, RMS force %.3f\n",
1321 gmx_step_str(mdstep
, sbuf
), number_steps
, df
[Min
]);
1324 /* Copy back the coordinates and the forces */
1325 std::copy(pos
[Min
].begin(), pos
[Min
].end(), state
->x
.begin());
1326 std::copy(force
[Min
].begin(), force
[Min
].end(), f
->begin());
1329 void done_shellfc(FILE *fplog
, gmx_shellfc_t
*shfc
, gmx_int64_t numSteps
)
1331 if (shfc
&& fplog
&& numSteps
> 0)
1333 double numStepsAsDouble
= static_cast<double>(numSteps
);
1334 fprintf(fplog
, "Fraction of iterations that converged: %.2f %%\n",
1335 (shfc
->numConvergedIterations
*100.0)/numStepsAsDouble
);
1336 fprintf(fplog
, "Average number of force evaluations per MD step: %.2f\n\n",
1337 shfc
->numForceEvaluations
/numStepsAsDouble
);
1340 // TODO Deallocate memory in shfc