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, 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/domdec.h"
50 #include "gromacs/domdec/domdec_struct.h"
51 #include "gromacs/gmxlib/chargegroup.h"
52 #include "gromacs/gmxlib/network.h"
53 #include "gromacs/math/functions.h"
54 #include "gromacs/math/units.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/math/vecdump.h"
57 #include "gromacs/mdlib/constr.h"
58 #include "gromacs/mdlib/force.h"
59 #include "gromacs/mdlib/mdrun.h"
60 #include "gromacs/mdlib/sim_util.h"
61 #include "gromacs/mdlib/vsite.h"
62 #include "gromacs/mdtypes/commrec.h"
63 #include "gromacs/mdtypes/inputrec.h"
64 #include "gromacs/mdtypes/md_enums.h"
65 #include "gromacs/pbcutil/mshift.h"
66 #include "gromacs/pbcutil/pbc.h"
67 #include "gromacs/topology/mtop_util.h"
68 #include "gromacs/utility/arraysize.h"
69 #include "gromacs/utility/cstringutil.h"
70 #include "gromacs/utility/fatalerror.h"
71 #include "gromacs/utility/smalloc.h"
75 int shell
; /* The shell id */
76 int nucl1
, nucl2
, nucl3
; /* The nuclei connected to the shell */
77 /* gmx_bool bInterCG; */ /* Coupled to nuclei outside cg? */
78 real k
; /* force constant */
79 real k_1
; /* 1 over force constant */
85 struct gmx_shellfc_t
{
86 int nshell_gl
; /* The number of shells in the system */
87 t_shell
*shell_gl
; /* All the shells (for DD only) */
88 int *shell_index_gl
; /* Global shell index (for DD only) */
89 gmx_bool bInterCG
; /* Are there inter charge-group shells? */
90 int nshell
; /* The number of local shells */
91 t_shell
*shell
; /* The local shells */
92 int shell_nalloc
; /* The allocation size of shell */
93 gmx_bool bPredict
; /* Predict shell positions */
94 gmx_bool bRequireInit
; /* Require initialization of shell positions */
95 int nflexcon
; /* The number of flexible constraints */
96 rvec
*x
[2]; /* Array for iterative minimization */
97 rvec
*f
[2]; /* Array for iterative minimization */
98 int x_nalloc
; /* The allocation size of x and f */
99 rvec
*acc_dir
; /* Acceleration direction for flexcon */
100 rvec
*x_old
; /* Old coordinates for flexcon */
101 int flex_nalloc
; /* The allocation size of acc_dir and x_old */
102 rvec
*adir_xnold
; /* Work space for init_adir */
103 rvec
*adir_xnew
; /* Work space for init_adir */
104 int adir_nalloc
; /* Work space for init_adir */
105 std::int64_t numForceEvaluations
; /* Total number of force evaluations */
106 int numConvergedIterations
; /* Total number of iterations that converged */
110 static void pr_shell(FILE *fplog
, int ns
, t_shell s
[])
114 fprintf(fplog
, "SHELL DATA\n");
115 fprintf(fplog
, "%5s %8s %5s %5s %5s\n",
116 "Shell", "Force k", "Nucl1", "Nucl2", "Nucl3");
117 for (i
= 0; (i
< ns
); i
++)
119 fprintf(fplog
, "%5d %8.3f %5d", s
[i
].shell
, 1.0/s
[i
].k_1
, s
[i
].nucl1
);
122 fprintf(fplog
, " %5d\n", s
[i
].nucl2
);
124 else if (s
[i
].nnucl
== 3)
126 fprintf(fplog
, " %5d %5d\n", s
[i
].nucl2
, s
[i
].nucl3
);
130 fprintf(fplog
, "\n");
135 /* TODO The remain call of this function passes non-NULL mass and NULL
136 * mtop, so this routine can be simplified.
138 * The other code path supported doing prediction before the MD loop
139 * started, but even when called, the prediction was always
140 * over-written by a subsequent call in the MD loop, so has been
142 static void predict_shells(FILE *fplog
, rvec x
[], rvec v
[], real dt
,
144 real mass
[], gmx_mtop_t
*mtop
, gmx_bool bInit
)
146 int i
, m
, s1
, n1
, n2
, n3
;
147 real dt_1
, fudge
, tm
, m1
, m2
, m3
;
149 gmx_mtop_atomlookup_t alook
= NULL
;
154 alook
= gmx_mtop_atomlookup_init(mtop
);
157 /* We introduce a fudge factor for performance reasons: with this choice
158 * the initial force on the shells is about a factor of two lower than
167 fprintf(fplog
, "RELAX: Using prediction for initial shell placement\n");
178 for (i
= 0; (i
< ns
); i
++)
189 for (m
= 0; (m
< DIM
); m
++)
191 x
[s1
][m
] += ptr
[n1
][m
]*dt_1
;
204 /* Not the correct masses with FE, but it is just a prediction... */
205 gmx_mtop_atomnr_to_atom(alook
, n1
, &atom
);
207 gmx_mtop_atomnr_to_atom(alook
, n2
, &atom
);
211 for (m
= 0; (m
< DIM
); m
++)
213 x
[s1
][m
] += (m1
*ptr
[n1
][m
]+m2
*ptr
[n2
][m
])*tm
;
228 /* Not the correct masses with FE, but it is just a prediction... */
229 gmx_mtop_atomnr_to_atom(alook
, n1
, &atom
);
231 gmx_mtop_atomnr_to_atom(alook
, n2
, &atom
);
233 gmx_mtop_atomnr_to_atom(alook
, n3
, &atom
);
236 tm
= dt_1
/(m1
+m2
+m3
);
237 for (m
= 0; (m
< DIM
); m
++)
239 x
[s1
][m
] += (m1
*ptr
[n1
][m
]+m2
*ptr
[n2
][m
]+m3
*ptr
[n3
][m
])*tm
;
243 gmx_fatal(FARGS
, "Shell %d has %d nuclei!", i
, s
[i
].nnucl
);
249 gmx_mtop_atomlookup_destroy(alook
);
253 /*! \brief Count the different particle types in a system
255 * Routine prints a warning to stderr in case an unknown particle type
257 * \param[in] fplog Print what we have found if not NULL
258 * \param[in] mtop Molecular topology.
259 * \returns Array holding the number of particles of a type
261 static std::array
<int, eptNR
> countPtypes(FILE *fplog
,
264 std::array
<int, eptNR
> nptype
= { { 0 } };
265 /* Count number of shells, and find their indices */
266 for (int i
= 0; (i
< eptNR
); i
++)
271 gmx_mtop_atomloop_block_t aloopb
= gmx_mtop_atomloop_block_init(mtop
);
274 while (gmx_mtop_atomloop_block_next(aloopb
, &atom
, &nmol
))
281 nptype
[atom
->ptype
] += nmol
;
284 fprintf(stderr
, "Warning unsupported particle type %d in countPtypes",
285 static_cast<int>(atom
->ptype
));
290 /* Print the number of each particle type */
292 for (const auto &i
: nptype
)
296 fprintf(fplog
, "There are: %d %ss\n", i
, ptype_str
[n
]);
304 gmx_shellfc_t
*init_shell_flexcon(FILE *fplog
,
305 gmx_mtop_t
*mtop
, int nflexcon
,
307 bool usingDomainDecomposition
)
311 int *shell_index
= NULL
, *at2cg
;
315 int i
, j
, type
, mb
, a_offset
, cg
, mol
, ftype
, nra
;
317 int aS
, aN
= 0; /* Shell and nucleus */
318 int bondtypes
[] = { F_BONDS
, F_HARMONIC
, F_CUBICBONDS
, F_POLARIZATION
, F_ANHARM_POL
, F_WATER_POL
};
319 #define NBT asize(bondtypes)
321 gmx_mtop_atomloop_all_t aloop
;
322 gmx_ffparams_t
*ffparams
;
323 gmx_molblock_t
*molb
;
327 std::array
<int, eptNR
> n
= countPtypes(fplog
, mtop
);
328 nshell
= n
[eptShell
];
330 if (nshell
== 0 && nflexcon
== 0)
332 /* We're not doing shells or flexible constraints */
337 shfc
->nflexcon
= nflexcon
;
339 if (nstcalcenergy
!= 1)
341 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
);
343 if (usingDomainDecomposition
)
345 gmx_fatal(FARGS
, "Shell particles are not implemented with domain decomposition, use a single rank");
353 /* We have shells: fill the shell data structure */
355 /* Global system sized array, this should be avoided */
356 snew(shell_index
, mtop
->natoms
);
358 aloop
= gmx_mtop_atomloop_all_init(mtop
);
360 while (gmx_mtop_atomloop_all_next(aloop
, &i
, &atom
))
362 if (atom
->ptype
== eptShell
)
364 shell_index
[i
] = nshell
++;
370 /* Initiate the shell structures */
371 for (i
= 0; (i
< nshell
); i
++)
378 /* shell[i].bInterCG=FALSE; */
383 ffparams
= &mtop
->ffparams
;
385 /* Now fill the structures */
386 shfc
->bInterCG
= FALSE
;
389 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
391 molb
= &mtop
->molblock
[mb
];
392 molt
= &mtop
->moltype
[molb
->type
];
395 snew(at2cg
, molt
->atoms
.nr
);
396 for (cg
= 0; cg
< cgs
->nr
; cg
++)
398 for (i
= cgs
->index
[cg
]; i
< cgs
->index
[cg
+1]; i
++)
404 atom
= molt
->atoms
.atom
;
405 for (mol
= 0; mol
< molb
->nmol
; mol
++)
407 for (j
= 0; (j
< NBT
); j
++)
409 ia
= molt
->ilist
[bondtypes
[j
]].iatoms
;
410 for (i
= 0; (i
< molt
->ilist
[bondtypes
[j
]].nr
); )
413 ftype
= ffparams
->functype
[type
];
414 nra
= interaction_function
[ftype
].nratoms
;
416 /* Check whether we have a bond with a shell */
419 switch (bondtypes
[j
])
426 if (atom
[ia
[1]].ptype
== eptShell
)
431 else if (atom
[ia
[2]].ptype
== eptShell
)
438 aN
= ia
[4]; /* Dummy */
439 aS
= ia
[5]; /* Shell */
442 gmx_fatal(FARGS
, "Death Horror: %s, %d", __FILE__
, __LINE__
);
449 /* Check whether one of the particles is a shell... */
450 nsi
= shell_index
[a_offset
+aS
];
451 if ((nsi
< 0) || (nsi
>= nshell
))
453 gmx_fatal(FARGS
, "nsi is %d should be within 0 - %d. aS = %d",
456 if (shell
[nsi
].shell
== -1)
458 shell
[nsi
].shell
= a_offset
+ aS
;
461 else if (shell
[nsi
].shell
!= a_offset
+aS
)
463 gmx_fatal(FARGS
, "Weird stuff in %s, %d", __FILE__
, __LINE__
);
466 if (shell
[nsi
].nucl1
== -1)
468 shell
[nsi
].nucl1
= a_offset
+ aN
;
470 else if (shell
[nsi
].nucl2
== -1)
472 shell
[nsi
].nucl2
= a_offset
+ aN
;
474 else if (shell
[nsi
].nucl3
== -1)
476 shell
[nsi
].nucl3
= a_offset
+ aN
;
482 pr_shell(fplog
, ns
, shell
);
484 gmx_fatal(FARGS
, "Can not handle more than three bonds per shell\n");
486 if (at2cg
[aS
] != at2cg
[aN
])
488 /* shell[nsi].bInterCG = TRUE; */
489 shfc
->bInterCG
= TRUE
;
492 switch (bondtypes
[j
])
496 shell
[nsi
].k
+= ffparams
->iparams
[type
].harmonic
.krA
;
499 shell
[nsi
].k
+= ffparams
->iparams
[type
].cubic
.kb
;
503 if (!gmx_within_tol(qS
, atom
[aS
].qB
, GMX_REAL_EPS
*10))
505 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);
507 shell
[nsi
].k
+= gmx::square(qS
)*ONE_4PI_EPS0
/
508 ffparams
->iparams
[type
].polarize
.alpha
;
511 if (!gmx_within_tol(qS
, atom
[aS
].qB
, GMX_REAL_EPS
*10))
513 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);
515 alpha
= (ffparams
->iparams
[type
].wpol
.al_x
+
516 ffparams
->iparams
[type
].wpol
.al_y
+
517 ffparams
->iparams
[type
].wpol
.al_z
)/3.0;
518 shell
[nsi
].k
+= gmx::square(qS
)*ONE_4PI_EPS0
/alpha
;
521 gmx_fatal(FARGS
, "Death Horror: %s, %d", __FILE__
, __LINE__
);
529 a_offset
+= molt
->atoms
.nr
;
531 /* Done with this molecule type */
535 /* Verify whether it's all correct */
538 gmx_fatal(FARGS
, "Something weird with shells. They may not be bonded to something");
541 for (i
= 0; (i
< ns
); i
++)
543 shell
[i
].k_1
= 1.0/shell
[i
].k
;
548 pr_shell(debug
, ns
, shell
);
552 shfc
->nshell_gl
= ns
;
553 shfc
->shell_gl
= shell
;
554 shfc
->shell_index_gl
= shell_index
;
556 shfc
->bPredict
= (getenv("GMX_NOPREDICT") == NULL
);
557 shfc
->bRequireInit
= FALSE
;
562 fprintf(fplog
, "\nWill never predict shell positions\n");
567 shfc
->bRequireInit
= (getenv("GMX_REQUIRE_SHELL_INIT") != NULL
);
568 if (shfc
->bRequireInit
&& fplog
)
570 fprintf(fplog
, "\nWill always initiate shell positions\n");
580 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");
582 /* Prediction improves performance, so we should implement either:
583 * 1. communication for the atoms needed for prediction
584 * 2. prediction using the velocities of shells; currently the
585 * shell velocities are zeroed, it's a bit tricky to keep
586 * track of the shell displacements and thus the velocity.
588 shfc
->bPredict
= FALSE
;
595 void make_local_shells(t_commrec
*cr
, t_mdatoms
*md
,
599 int a0
, a1
, *ind
, nshell
, i
;
600 gmx_domdec_t
*dd
= NULL
;
602 if (DOMAINDECOMP(cr
))
610 /* Single node: we need all shells, just copy the pointer */
611 shfc
->nshell
= shfc
->nshell_gl
;
612 shfc
->shell
= shfc
->shell_gl
;
617 ind
= shfc
->shell_index_gl
;
621 for (i
= a0
; i
< a1
; i
++)
623 if (md
->ptype
[i
] == eptShell
)
625 if (nshell
+1 > shfc
->shell_nalloc
)
627 shfc
->shell_nalloc
= over_alloc_dd(nshell
+1);
628 srenew(shell
, shfc
->shell_nalloc
);
632 shell
[nshell
] = shfc
->shell_gl
[ind
[dd
->gatindex
[i
]]];
636 shell
[nshell
] = shfc
->shell_gl
[ind
[i
]];
639 /* With inter-cg shells we can no do shell prediction,
640 * so we do not need the nuclei numbers.
644 shell
[nshell
].nucl1
= i
+ shell
[nshell
].nucl1
- shell
[nshell
].shell
;
645 if (shell
[nshell
].nnucl
> 1)
647 shell
[nshell
].nucl2
= i
+ shell
[nshell
].nucl2
- shell
[nshell
].shell
;
649 if (shell
[nshell
].nnucl
> 2)
651 shell
[nshell
].nucl3
= i
+ shell
[nshell
].nucl3
- shell
[nshell
].shell
;
654 shell
[nshell
].shell
= i
;
659 shfc
->nshell
= nshell
;
663 static void do_1pos(rvec xnew
, rvec xold
, rvec f
, real step
)
681 static void do_1pos3(rvec xnew
, rvec xold
, rvec f
, rvec step
)
699 static void directional_sd(rvec xold
[], rvec xnew
[], rvec acc_dir
[],
700 int start
, int homenr
, real step
)
704 for (i
= start
; i
< homenr
; i
++)
706 do_1pos(xnew
[i
], xold
[i
], acc_dir
[i
], step
);
710 static void shell_pos_sd(rvec xcur
[], rvec xnew
[], rvec f
[],
711 int ns
, t_shell s
[], int count
)
713 const real step_scale_min
= 0.8,
714 step_scale_increment
= 0.2,
715 step_scale_max
= 1.2,
716 step_scale_multiple
= (step_scale_max
- step_scale_min
) / step_scale_increment
;
721 real step_min
, step_max
;
726 for (i
= 0; (i
< ns
); i
++)
731 for (d
= 0; d
< DIM
; d
++)
733 s
[i
].step
[d
] = s
[i
].k_1
;
735 step_min
= std::min(step_min
, s
[i
].step
[d
]);
736 step_max
= std::max(step_max
, s
[i
].step
[d
]);
742 for (d
= 0; d
< DIM
; d
++)
744 dx
= xcur
[shell
][d
] - s
[i
].xold
[d
];
745 df
= f
[shell
][d
] - s
[i
].fold
[d
];
746 /* -dx/df gets used to generate an interpolated value, but would
747 * cause a NaN if df were binary-equal to zero. Values close to
748 * zero won't cause problems (because of the min() and max()), so
749 * just testing for binary inequality is OK. */
753 /* Scale the step size by a factor interpolated from
754 * step_scale_min to step_scale_max, as k_est goes from 0 to
755 * step_scale_multiple * s[i].step[d] */
757 step_scale_min
* s
[i
].step
[d
] +
758 step_scale_increment
* std::min(step_scale_multiple
* s
[i
].step
[d
], std::max(k_est
, zero
));
763 if (gmx_numzero(dx
)) /* 0 == dx */
765 /* Likely this will never happen, but if it does just
766 * don't scale the step. */
770 s
[i
].step
[d
] *= step_scale_max
;
774 step_min
= std::min(step_min
, s
[i
].step
[d
]);
775 step_max
= std::max(step_max
, s
[i
].step
[d
]);
779 copy_rvec(xcur
[shell
], s
[i
].xold
);
780 copy_rvec(f
[shell
], s
[i
].fold
);
782 do_1pos3(xnew
[shell
], xcur
[shell
], f
[shell
], s
[i
].step
);
786 fprintf(debug
, "shell[%d] = %d\n", i
, shell
);
787 pr_rvec(debug
, 0, "fshell", f
[shell
], DIM
, TRUE
);
788 pr_rvec(debug
, 0, "xold", xcur
[shell
], DIM
, TRUE
);
789 pr_rvec(debug
, 0, "step", s
[i
].step
, DIM
, TRUE
);
790 pr_rvec(debug
, 0, "xnew", xnew
[shell
], DIM
, TRUE
);
794 printf("step %.3e %.3e\n", step_min
, step_max
);
798 static void decrease_step_size(int nshell
, t_shell s
[])
802 for (i
= 0; i
< nshell
; i
++)
804 svmul(0.8, s
[i
].step
, s
[i
].step
);
808 static void print_epot(FILE *fp
, gmx_int64_t mdstep
, int count
, real epot
, real df
,
809 int ndir
, real sf_dir
)
813 fprintf(fp
, "MDStep=%5s/%2d EPot: %12.8e, rmsF: %6.2e",
814 gmx_step_str(mdstep
, buf
), count
, epot
, df
);
817 fprintf(fp
, ", dir. rmsF: %6.2e\n", std::sqrt(sf_dir
/ndir
));
826 static real
rms_force(t_commrec
*cr
, rvec f
[], int ns
, t_shell s
[],
827 int ndir
, real
*sf_dir
, real
*Epot
)
833 for (i
= 0; i
< ns
; i
++)
836 buf
[0] += norm2(f
[shell
]);
845 gmx_sumd(4, buf
, cr
);
846 ntot
= (int)(buf
[1] + 0.5);
852 return (ntot
? std::sqrt(buf
[0]/ntot
) : 0);
855 static void check_pbc(FILE *fp
, rvec x
[], int shell
)
860 for (m
= 0; (m
< DIM
); m
++)
862 if (fabs(x
[shell
][m
]-x
[now
][m
]) > 0.3)
864 pr_rvecs(fp
, 0, "SHELL-X", x
+now
, 5);
870 static void dump_shells(FILE *fp
, rvec x
[], rvec f
[], real ftol
, int ns
, t_shell s
[])
875 ft2
= gmx::square(ftol
);
877 for (i
= 0; (i
< ns
); i
++)
880 ff2
= iprod(f
[shell
], f
[shell
]);
883 fprintf(fp
, "SHELL %5d, force %10.5f %10.5f %10.5f, |f| %10.5f\n",
884 shell
, f
[shell
][XX
], f
[shell
][YY
], f
[shell
][ZZ
], std::sqrt(ff2
));
886 check_pbc(fp
, x
, shell
);
890 static void init_adir(FILE *log
, gmx_shellfc_t
*shfc
,
891 gmx_constr_t constr
, t_idef
*idef
, t_inputrec
*ir
,
892 t_commrec
*cr
, int dd_ac1
,
893 gmx_int64_t step
, t_mdatoms
*md
, int start
, int end
,
894 rvec
*x_old
, rvec
*x_init
, rvec
*x
,
895 rvec
*f
, rvec
*acc_dir
,
896 gmx_bool bMolPBC
, matrix box
,
897 real
*lambda
, real
*dvdlambda
, t_nrnb
*nrnb
)
902 unsigned short *ptype
;
904 if (DOMAINDECOMP(cr
))
912 if (n
> shfc
->adir_nalloc
)
914 shfc
->adir_nalloc
= over_alloc_dd(n
);
915 srenew(shfc
->adir_xnold
, shfc
->adir_nalloc
);
916 srenew(shfc
->adir_xnew
, shfc
->adir_nalloc
);
918 xnold
= shfc
->adir_xnold
;
919 xnew
= shfc
->adir_xnew
;
925 /* Does NOT work with freeze or acceleration groups (yet) */
926 for (n
= start
; n
< end
; n
++)
928 w_dt
= md
->invmass
[n
]*dt
;
930 for (d
= 0; d
< DIM
; d
++)
932 if ((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
))
934 xnold
[n
-start
][d
] = x
[n
][d
] - (x_init
[n
][d
] - x_old
[n
][d
]);
935 xnew
[n
-start
][d
] = 2*x
[n
][d
] - x_old
[n
][d
] + f
[n
][d
]*w_dt
*dt
;
939 xnold
[n
-start
][d
] = x
[n
][d
];
940 xnew
[n
-start
][d
] = x
[n
][d
];
944 constrain(log
, FALSE
, FALSE
, constr
, idef
, ir
, cr
, step
, 0, 1.0, md
,
945 x
, xnold
-start
, NULL
, bMolPBC
, box
,
946 lambda
[efptBONDED
], &(dvdlambda
[efptBONDED
]),
947 NULL
, NULL
, nrnb
, econqCoord
);
948 constrain(log
, FALSE
, FALSE
, constr
, idef
, ir
, cr
, step
, 0, 1.0, md
,
949 x
, xnew
-start
, NULL
, bMolPBC
, box
,
950 lambda
[efptBONDED
], &(dvdlambda
[efptBONDED
]),
951 NULL
, NULL
, nrnb
, econqCoord
);
953 for (n
= start
; n
< end
; n
++)
955 for (d
= 0; d
< DIM
; d
++)
958 -(2*x
[n
][d
]-xnold
[n
-start
][d
]-xnew
[n
-start
][d
])/gmx::square(dt
)
959 - f
[n
][d
]*md
->invmass
[n
];
961 clear_rvec(acc_dir
[n
]);
964 /* Project the acceleration on the old bond directions */
965 constrain(log
, FALSE
, FALSE
, constr
, idef
, ir
, cr
, step
, 0, 1.0, md
,
966 x_old
, xnew
-start
, acc_dir
, bMolPBC
, box
,
967 lambda
[efptBONDED
], &(dvdlambda
[efptBONDED
]),
968 NULL
, NULL
, nrnb
, econqDeriv_FlexCon
);
971 void relax_shell_flexcon(FILE *fplog
, t_commrec
*cr
, gmx_bool bVerbose
,
972 gmx_int64_t mdstep
, t_inputrec
*inputrec
,
973 gmx_bool bDoNS
, int force_flags
,
976 gmx_enerdata_t
*enerd
, t_fcdata
*fcd
,
977 t_state
*state
, rvec f
[],
980 t_nrnb
*nrnb
, gmx_wallcycle_t wcycle
,
982 gmx_groups_t
*groups
,
986 double t
, rvec mu_tot
,
993 rvec
*pos
[2], *force
[2], *acc_dir
= NULL
, *x_old
= NULL
;
998 gmx_bool bCont
, bInit
, bConverged
;
999 int nat
, dd_ac0
, dd_ac1
= 0, i
;
1000 int start
= 0, homenr
= md
->homenr
, end
= start
+homenr
, cg0
, cg1
;
1001 int nflexcon
, number_steps
, d
, Min
= 0, count
= 0;
1002 #define Try (1-Min) /* At start Try = 1 */
1004 bCont
= (mdstep
== inputrec
->init_step
) && inputrec
->bContinuation
;
1005 bInit
= (mdstep
== inputrec
->init_step
) || shfc
->bRequireInit
;
1006 ftol
= inputrec
->em_tol
;
1007 number_steps
= inputrec
->niter
;
1008 nshell
= shfc
->nshell
;
1009 shell
= shfc
->shell
;
1010 nflexcon
= shfc
->nflexcon
;
1014 if (DOMAINDECOMP(cr
))
1016 nat
= dd_natoms_vsite(cr
->dd
);
1019 dd_get_constraint_range(cr
->dd
, &dd_ac0
, &dd_ac1
);
1020 nat
= std::max(nat
, dd_ac1
);
1025 nat
= state
->natoms
;
1028 if (nat
> shfc
->x_nalloc
)
1030 /* Allocate local arrays */
1031 shfc
->x_nalloc
= over_alloc_dd(nat
);
1032 for (i
= 0; (i
< 2); i
++)
1034 srenew(shfc
->x
[i
], shfc
->x_nalloc
);
1035 srenew(shfc
->f
[i
], shfc
->x_nalloc
);
1038 for (i
= 0; (i
< 2); i
++)
1040 pos
[i
] = shfc
->x
[i
];
1041 force
[i
] = shfc
->f
[i
];
1044 if (bDoNS
&& inputrec
->ePBC
!= epbcNONE
&& !DOMAINDECOMP(cr
))
1046 /* This is the only time where the coordinates are used
1047 * before do_force is called, which normally puts all
1048 * charge groups in the box.
1050 if (inputrec
->cutoff_scheme
== ecutsVERLET
)
1052 put_atoms_in_box_omp(fr
->ePBC
, state
->box
, md
->homenr
, state
->x
);
1058 put_charge_groups_in_box(fplog
, cg0
, cg1
, fr
->ePBC
, state
->box
,
1059 &(top
->cgs
), state
->x
, fr
->cg_cm
);
1064 mk_mshift(fplog
, graph
, fr
->ePBC
, state
->box
, state
->x
);
1068 /* After this all coordinate arrays will contain whole charge groups */
1071 shift_self(graph
, state
->box
, state
->x
);
1076 if (nat
> shfc
->flex_nalloc
)
1078 shfc
->flex_nalloc
= over_alloc_dd(nat
);
1079 srenew(shfc
->acc_dir
, shfc
->flex_nalloc
);
1080 srenew(shfc
->x_old
, shfc
->flex_nalloc
);
1082 acc_dir
= shfc
->acc_dir
;
1083 x_old
= shfc
->x_old
;
1084 for (i
= 0; i
< homenr
; i
++)
1086 for (d
= 0; d
< DIM
; d
++)
1089 state
->x
[start
+i
][d
] - state
->v
[start
+i
][d
]*inputrec
->delta_t
;
1094 /* Do a prediction of the shell positions */
1095 if (shfc
->bPredict
&& !bCont
)
1097 predict_shells(fplog
, state
->x
, state
->v
, inputrec
->delta_t
, nshell
, shell
,
1098 md
->massT
, NULL
, bInit
);
1101 /* do_force expected the charge groups to be in the box */
1104 unshift_self(graph
, state
->box
, state
->x
);
1107 /* Calculate the forces first time around */
1110 pr_rvecs(debug
, 0, "x b4 do_force", state
->x
+ start
, homenr
);
1112 do_force(fplog
, cr
, inputrec
, mdstep
, nrnb
, wcycle
, top
, groups
,
1113 state
->box
, state
->x
, &state
->hist
,
1114 force
[Min
], force_vir
, md
, enerd
, fcd
,
1115 state
->lambda
, graph
,
1116 fr
, vsite
, mu_tot
, t
, fp_field
, NULL
, bBornRadii
,
1117 (bDoNS
? GMX_FORCE_NS
: 0) | force_flags
);
1122 init_adir(fplog
, shfc
,
1123 constr
, idef
, inputrec
, cr
, dd_ac1
, mdstep
, md
, start
, end
,
1124 shfc
->x_old
-start
, state
->x
, state
->x
, force
[Min
],
1125 shfc
->acc_dir
-start
,
1126 fr
->bMolPBC
, state
->box
, state
->lambda
, &dum
, nrnb
);
1128 for (i
= start
; i
< end
; i
++)
1130 sf_dir
+= md
->massT
[i
]*norm2(shfc
->acc_dir
[i
-start
]);
1134 Epot
[Min
] = enerd
->term
[F_EPOT
];
1136 df
[Min
] = rms_force(cr
, shfc
->f
[Min
], nshell
, shell
, nflexcon
, &sf_dir
, &Epot
[Min
]);
1140 fprintf(debug
, "df = %g %g\n", df
[Min
], df
[Try
]);
1145 pr_rvecs(debug
, 0, "force0", force
[Min
], md
->nr
);
1148 if (nshell
+nflexcon
> 0)
1150 /* Copy x to pos[Min] & pos[Try]: during minimization only the
1151 * shell positions are updated, therefore the other particles must
1154 memcpy(pos
[Min
], state
->x
, nat
*sizeof(state
->x
[0]));
1155 memcpy(pos
[Try
], state
->x
, nat
*sizeof(state
->x
[0]));
1158 if (bVerbose
&& MASTER(cr
))
1160 print_epot(stdout
, mdstep
, 0, Epot
[Min
], df
[Min
], nflexcon
, sf_dir
);
1165 fprintf(debug
, "%17s: %14.10e\n",
1166 interaction_function
[F_EKIN
].longname
, enerd
->term
[F_EKIN
]);
1167 fprintf(debug
, "%17s: %14.10e\n",
1168 interaction_function
[F_EPOT
].longname
, enerd
->term
[F_EPOT
]);
1169 fprintf(debug
, "%17s: %14.10e\n",
1170 interaction_function
[F_ETOT
].longname
, enerd
->term
[F_ETOT
]);
1171 fprintf(debug
, "SHELLSTEP %s\n", gmx_step_str(mdstep
, sbuf
));
1174 /* First check whether we should do shells, or whether the force is
1175 * low enough even without minimization.
1177 bConverged
= (df
[Min
] < ftol
);
1179 for (count
= 1; (!(bConverged
) && (count
< number_steps
)); count
++)
1183 construct_vsites(vsite
, pos
[Min
], inputrec
->delta_t
, state
->v
,
1184 idef
->iparams
, idef
->il
,
1185 fr
->ePBC
, fr
->bMolPBC
, cr
, state
->box
);
1190 init_adir(fplog
, shfc
,
1191 constr
, idef
, inputrec
, cr
, dd_ac1
, mdstep
, md
, start
, end
,
1192 x_old
-start
, state
->x
, pos
[Min
], force
[Min
], acc_dir
-start
,
1193 fr
->bMolPBC
, state
->box
, state
->lambda
, &dum
, nrnb
);
1195 directional_sd(pos
[Min
], pos
[Try
], acc_dir
-start
, start
, end
,
1199 /* New positions, Steepest descent */
1200 shell_pos_sd(pos
[Min
], pos
[Try
], force
[Min
], nshell
, shell
, count
);
1202 /* do_force expected the charge groups to be in the box */
1205 unshift_self(graph
, state
->box
, pos
[Try
]);
1210 pr_rvecs(debug
, 0, "RELAX: pos[Min] ", pos
[Min
] + start
, homenr
);
1211 pr_rvecs(debug
, 0, "RELAX: pos[Try] ", pos
[Try
] + start
, homenr
);
1213 /* Try the new positions */
1214 do_force(fplog
, cr
, inputrec
, 1, nrnb
, wcycle
,
1215 top
, groups
, state
->box
, pos
[Try
], &state
->hist
,
1216 force
[Try
], force_vir
,
1217 md
, enerd
, fcd
, state
->lambda
, graph
,
1218 fr
, vsite
, mu_tot
, t
, fp_field
, NULL
, bBornRadii
,
1223 pr_rvecs(debug
, 0, "RELAX: force[Min]", force
[Min
] + start
, homenr
);
1224 pr_rvecs(debug
, 0, "RELAX: force[Try]", force
[Try
] + start
, homenr
);
1229 init_adir(fplog
, shfc
,
1230 constr
, idef
, inputrec
, cr
, dd_ac1
, mdstep
, md
, start
, end
,
1231 x_old
-start
, state
->x
, pos
[Try
], force
[Try
], acc_dir
-start
,
1232 fr
->bMolPBC
, state
->box
, state
->lambda
, &dum
, nrnb
);
1234 for (i
= start
; i
< end
; i
++)
1236 sf_dir
+= md
->massT
[i
]*norm2(acc_dir
[i
-start
]);
1240 Epot
[Try
] = enerd
->term
[F_EPOT
];
1242 df
[Try
] = rms_force(cr
, force
[Try
], nshell
, shell
, nflexcon
, &sf_dir
, &Epot
[Try
]);
1246 fprintf(debug
, "df = %g %g\n", df
[Min
], df
[Try
]);
1253 pr_rvecs(debug
, 0, "F na do_force", force
[Try
] + start
, homenr
);
1257 fprintf(debug
, "SHELL ITER %d\n", count
);
1258 dump_shells(debug
, pos
[Try
], force
[Try
], ftol
, nshell
, shell
);
1262 if (bVerbose
&& MASTER(cr
))
1264 print_epot(stdout
, mdstep
, count
, Epot
[Try
], df
[Try
], nflexcon
, sf_dir
);
1267 bConverged
= (df
[Try
] < ftol
);
1269 if ((df
[Try
] < df
[Min
]))
1273 fprintf(debug
, "Swapping Min and Try\n");
1277 /* Correct the velocities for the flexible constraints */
1278 invdt
= 1/inputrec
->delta_t
;
1279 for (i
= start
; i
< end
; i
++)
1281 for (d
= 0; d
< DIM
; d
++)
1283 state
->v
[i
][d
] += (pos
[Try
][i
][d
] - pos
[Min
][i
][d
])*invdt
;
1291 decrease_step_size(nshell
, shell
);
1294 shfc
->numForceEvaluations
+= count
;
1297 shfc
->numConvergedIterations
++;
1299 if (MASTER(cr
) && !(bConverged
))
1301 /* Note that the energies and virial are incorrect when not converged */
1305 "step %s: EM did not converge in %d iterations, RMS force %.3f\n",
1306 gmx_step_str(mdstep
, sbuf
), number_steps
, df
[Min
]);
1309 "step %s: EM did not converge in %d iterations, RMS force %.3f\n",
1310 gmx_step_str(mdstep
, sbuf
), number_steps
, df
[Min
]);
1313 /* Copy back the coordinates and the forces */
1314 memcpy(state
->x
, pos
[Min
], nat
*sizeof(state
->x
[0]));
1315 memcpy(f
, force
[Min
], nat
*sizeof(f
[0]));
1318 void done_shellfc(FILE *fplog
, gmx_shellfc_t
*shfc
, gmx_int64_t numSteps
)
1320 if (shfc
&& fplog
&& numSteps
> 0)
1322 double numStepsAsDouble
= static_cast<double>(numSteps
);
1323 fprintf(fplog
, "Fraction of iterations that converged: %.2f %%\n",
1324 (shfc
->numConvergedIterations
*100.0)/numStepsAsDouble
);
1325 fprintf(fplog
, "Average number of force evaluations per MD step: %.2f\n\n",
1326 shfc
->numForceEvaluations
/numStepsAsDouble
);
1329 // TODO Deallocate memory in shfc