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, 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.
43 #include "types/commrec.h"
44 #include "gromacs/utility/smalloc.h"
45 #include "gmx_fatal.h"
58 #include "mtop_util.h"
59 #include "chargegroup.h"
65 atom_id shell
; /* The shell id */
66 atom_id nucl1
, nucl2
, nucl3
; /* The nuclei connected to the shell */
67 /* gmx_bool bInterCG; */ /* Coupled to nuclei outside cg? */
68 real k
; /* force constant */
69 real k_1
; /* 1 over force constant */
75 typedef struct gmx_shellfc
{
76 int nshell_gl
; /* The number of shells in the system */
77 t_shell
*shell_gl
; /* All the shells (for DD only) */
78 int *shell_index_gl
; /* Global shell index (for DD only) */
79 gmx_bool bInterCG
; /* Are there inter charge-group shells? */
80 int nshell
; /* The number of local shells */
81 t_shell
*shell
; /* The local shells */
82 int shell_nalloc
; /* The allocation size of shell */
83 gmx_bool bPredict
; /* Predict shell positions */
84 gmx_bool bRequireInit
; /* Require initialization of shell positions */
85 int nflexcon
; /* The number of flexible constraints */
86 rvec
*x
[2]; /* Array for iterative minimization */
87 rvec
*f
[2]; /* Array for iterative minimization */
88 int x_nalloc
; /* The allocation size of x and f */
89 rvec
*acc_dir
; /* Acceleration direction for flexcon */
90 rvec
*x_old
; /* Old coordinates for flexcon */
91 int flex_nalloc
; /* The allocation size of acc_dir and x_old */
92 rvec
*adir_xnold
; /* Work space for init_adir */
93 rvec
*adir_xnew
; /* Work space for init_adir */
94 int adir_nalloc
; /* Work space for init_adir */
98 static void pr_shell(FILE *fplog
, int ns
, t_shell s
[])
102 fprintf(fplog
, "SHELL DATA\n");
103 fprintf(fplog
, "%5s %8s %5s %5s %5s\n",
104 "Shell", "Force k", "Nucl1", "Nucl2", "Nucl3");
105 for (i
= 0; (i
< ns
); i
++)
107 fprintf(fplog
, "%5d %8.3f %5d", s
[i
].shell
, 1.0/s
[i
].k_1
, s
[i
].nucl1
);
110 fprintf(fplog
, " %5d\n", s
[i
].nucl2
);
112 else if (s
[i
].nnucl
== 3)
114 fprintf(fplog
, " %5d %5d\n", s
[i
].nucl2
, s
[i
].nucl3
);
118 fprintf(fplog
, "\n");
123 static void predict_shells(FILE *fplog
, rvec x
[], rvec v
[], real dt
,
125 real mass
[], gmx_mtop_t
*mtop
, gmx_bool bInit
)
127 int i
, m
, s1
, n1
, n2
, n3
;
128 real dt_1
, dt_2
, dt_3
, fudge
, tm
, m1
, m2
, m3
;
130 gmx_mtop_atomlookup_t alook
= NULL
;
135 alook
= gmx_mtop_atomlookup_init(mtop
);
138 /* We introduce a fudge factor for performance reasons: with this choice
139 * the initial force on the shells is about a factor of two lower than
148 fprintf(fplog
, "RELAX: Using prediction for initial shell placement\n");
159 for (i
= 0; (i
< ns
); i
++)
170 for (m
= 0; (m
< DIM
); m
++)
172 x
[s1
][m
] += ptr
[n1
][m
]*dt_1
;
185 /* Not the correct masses with FE, but it is just a prediction... */
190 for (m
= 0; (m
< DIM
); m
++)
192 x
[s1
][m
] += (m1
*ptr
[n1
][m
]+m2
*ptr
[n2
][m
])*tm
;
207 /* Not the correct masses with FE, but it is just a prediction... */
208 gmx_mtop_atomnr_to_atom(alook
, n1
, &atom
);
210 gmx_mtop_atomnr_to_atom(alook
, n2
, &atom
);
212 gmx_mtop_atomnr_to_atom(alook
, n3
, &atom
);
215 tm
= dt_1
/(m1
+m2
+m3
);
216 for (m
= 0; (m
< DIM
); m
++)
218 x
[s1
][m
] += (m1
*ptr
[n1
][m
]+m2
*ptr
[n2
][m
]+m3
*ptr
[n3
][m
])*tm
;
222 gmx_fatal(FARGS
, "Shell %d has %d nuclei!", i
, s
[i
].nnucl
);
228 gmx_mtop_atomlookup_destroy(alook
);
232 gmx_shellfc_t
init_shell_flexcon(FILE *fplog
,
233 gmx_mtop_t
*mtop
, int nflexcon
,
236 struct gmx_shellfc
*shfc
;
238 int *shell_index
= NULL
, *at2cg
;
240 int n
[eptNR
], ns
, nshell
, nsi
;
241 int i
, j
, nmol
, type
, mb
, mt
, a_offset
, cg
, mol
, ftype
, nra
;
243 int aS
, aN
= 0; /* Shell and nucleus */
244 int bondtypes
[] = { F_BONDS
, F_HARMONIC
, F_CUBICBONDS
, F_POLARIZATION
, F_ANHARM_POL
, F_WATER_POL
};
245 #define NBT asize(bondtypes)
247 gmx_mtop_atomloop_block_t aloopb
;
248 gmx_mtop_atomloop_all_t aloop
;
249 gmx_ffparams_t
*ffparams
;
250 gmx_molblock_t
*molb
;
254 /* Count number of shells, and find their indices */
255 for (i
= 0; (i
< eptNR
); i
++)
260 aloopb
= gmx_mtop_atomloop_block_init(mtop
);
261 while (gmx_mtop_atomloop_block_next(aloopb
, &atom
, &nmol
))
263 n
[atom
->ptype
] += nmol
;
268 /* Print the number of each particle type */
269 for (i
= 0; (i
< eptNR
); i
++)
273 fprintf(fplog
, "There are: %d %ss\n", n
[i
], ptype_str
[i
]);
278 nshell
= n
[eptShell
];
280 if (nshell
== 0 && nflexcon
== 0)
282 /* We're not doing shells or flexible constraints */
287 shfc
->nflexcon
= nflexcon
;
294 /* We have shells: fill the shell data structure */
296 /* Global system sized array, this should be avoided */
297 snew(shell_index
, mtop
->natoms
);
299 aloop
= gmx_mtop_atomloop_all_init(mtop
);
301 while (gmx_mtop_atomloop_all_next(aloop
, &i
, &atom
))
303 if (atom
->ptype
== eptShell
)
305 shell_index
[i
] = nshell
++;
311 /* Initiate the shell structures */
312 for (i
= 0; (i
< nshell
); i
++)
314 shell
[i
].shell
= NO_ATID
;
316 shell
[i
].nucl1
= NO_ATID
;
317 shell
[i
].nucl2
= NO_ATID
;
318 shell
[i
].nucl3
= NO_ATID
;
319 /* shell[i].bInterCG=FALSE; */
324 ffparams
= &mtop
->ffparams
;
326 /* Now fill the structures */
327 shfc
->bInterCG
= FALSE
;
330 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
332 molb
= &mtop
->molblock
[mb
];
333 molt
= &mtop
->moltype
[molb
->type
];
336 snew(at2cg
, molt
->atoms
.nr
);
337 for (cg
= 0; cg
< cgs
->nr
; cg
++)
339 for (i
= cgs
->index
[cg
]; i
< cgs
->index
[cg
+1]; i
++)
345 atom
= molt
->atoms
.atom
;
346 for (mol
= 0; mol
< molb
->nmol
; mol
++)
348 for (j
= 0; (j
< NBT
); j
++)
350 ia
= molt
->ilist
[bondtypes
[j
]].iatoms
;
351 for (i
= 0; (i
< molt
->ilist
[bondtypes
[j
]].nr
); )
354 ftype
= ffparams
->functype
[type
];
355 nra
= interaction_function
[ftype
].nratoms
;
357 /* Check whether we have a bond with a shell */
360 switch (bondtypes
[j
])
367 if (atom
[ia
[1]].ptype
== eptShell
)
372 else if (atom
[ia
[2]].ptype
== eptShell
)
379 aN
= ia
[4]; /* Dummy */
380 aS
= ia
[5]; /* Shell */
383 gmx_fatal(FARGS
, "Death Horror: %s, %d", __FILE__
, __LINE__
);
390 /* Check whether one of the particles is a shell... */
391 nsi
= shell_index
[a_offset
+aS
];
392 if ((nsi
< 0) || (nsi
>= nshell
))
394 gmx_fatal(FARGS
, "nsi is %d should be within 0 - %d. aS = %d",
397 if (shell
[nsi
].shell
== NO_ATID
)
399 shell
[nsi
].shell
= a_offset
+ aS
;
402 else if (shell
[nsi
].shell
!= a_offset
+aS
)
404 gmx_fatal(FARGS
, "Weird stuff in %s, %d", __FILE__
, __LINE__
);
407 if (shell
[nsi
].nucl1
== NO_ATID
)
409 shell
[nsi
].nucl1
= a_offset
+ aN
;
411 else if (shell
[nsi
].nucl2
== NO_ATID
)
413 shell
[nsi
].nucl2
= a_offset
+ aN
;
415 else if (shell
[nsi
].nucl3
== NO_ATID
)
417 shell
[nsi
].nucl3
= a_offset
+ aN
;
423 pr_shell(fplog
, ns
, shell
);
425 gmx_fatal(FARGS
, "Can not handle more than three bonds per shell\n");
427 if (at2cg
[aS
] != at2cg
[aN
])
429 /* shell[nsi].bInterCG = TRUE; */
430 shfc
->bInterCG
= TRUE
;
433 switch (bondtypes
[j
])
437 shell
[nsi
].k
+= ffparams
->iparams
[type
].harmonic
.krA
;
440 shell
[nsi
].k
+= ffparams
->iparams
[type
].cubic
.kb
;
444 if (!gmx_within_tol(qS
, atom
[aS
].qB
, GMX_REAL_EPS
*10))
446 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);
448 shell
[nsi
].k
+= sqr(qS
)*ONE_4PI_EPS0
/
449 ffparams
->iparams
[type
].polarize
.alpha
;
452 if (!gmx_within_tol(qS
, atom
[aS
].qB
, GMX_REAL_EPS
*10))
454 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);
456 alpha
= (ffparams
->iparams
[type
].wpol
.al_x
+
457 ffparams
->iparams
[type
].wpol
.al_y
+
458 ffparams
->iparams
[type
].wpol
.al_z
)/3.0;
459 shell
[nsi
].k
+= sqr(qS
)*ONE_4PI_EPS0
/alpha
;
462 gmx_fatal(FARGS
, "Death Horror: %s, %d", __FILE__
, __LINE__
);
470 a_offset
+= molt
->atoms
.nr
;
472 /* Done with this molecule type */
476 /* Verify whether it's all correct */
479 gmx_fatal(FARGS
, "Something weird with shells. They may not be bonded to something");
482 for (i
= 0; (i
< ns
); i
++)
484 shell
[i
].k_1
= 1.0/shell
[i
].k
;
489 pr_shell(debug
, ns
, shell
);
493 shfc
->nshell_gl
= ns
;
494 shfc
->shell_gl
= shell
;
495 shfc
->shell_index_gl
= shell_index
;
497 shfc
->bPredict
= (getenv("GMX_NOPREDICT") == NULL
);
498 shfc
->bRequireInit
= FALSE
;
503 fprintf(fplog
, "\nWill never predict shell positions\n");
508 shfc
->bRequireInit
= (getenv("GMX_REQUIRE_SHELL_INIT") != NULL
);
509 if (shfc
->bRequireInit
&& fplog
)
511 fprintf(fplog
, "\nWill always initiate shell positions\n");
519 predict_shells(fplog
, x
, NULL
, 0, shfc
->nshell_gl
, shfc
->shell_gl
,
527 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");
529 /* Prediction improves performance, so we should implement either:
530 * 1. communication for the atoms needed for prediction
531 * 2. prediction using the velocities of shells; currently the
532 * shell velocities are zeroed, it's a bit tricky to keep
533 * track of the shell displacements and thus the velocity.
535 shfc
->bPredict
= FALSE
;
542 void make_local_shells(t_commrec
*cr
, t_mdatoms
*md
,
543 struct gmx_shellfc
*shfc
)
546 int a0
, a1
, *ind
, nshell
, i
;
547 gmx_domdec_t
*dd
= NULL
;
549 if (DOMAINDECOMP(cr
))
557 /* Single node: we need all shells, just copy the pointer */
558 shfc
->nshell
= shfc
->nshell_gl
;
559 shfc
->shell
= shfc
->shell_gl
;
564 ind
= shfc
->shell_index_gl
;
568 for (i
= a0
; i
< a1
; i
++)
570 if (md
->ptype
[i
] == eptShell
)
572 if (nshell
+1 > shfc
->shell_nalloc
)
574 shfc
->shell_nalloc
= over_alloc_dd(nshell
+1);
575 srenew(shell
, shfc
->shell_nalloc
);
579 shell
[nshell
] = shfc
->shell_gl
[ind
[dd
->gatindex
[i
]]];
583 shell
[nshell
] = shfc
->shell_gl
[ind
[i
]];
586 /* With inter-cg shells we can no do shell prediction,
587 * so we do not need the nuclei numbers.
591 shell
[nshell
].nucl1
= i
+ shell
[nshell
].nucl1
- shell
[nshell
].shell
;
592 if (shell
[nshell
].nnucl
> 1)
594 shell
[nshell
].nucl2
= i
+ shell
[nshell
].nucl2
- shell
[nshell
].shell
;
596 if (shell
[nshell
].nnucl
> 2)
598 shell
[nshell
].nucl3
= i
+ shell
[nshell
].nucl3
- shell
[nshell
].shell
;
601 shell
[nshell
].shell
= i
;
606 shfc
->nshell
= nshell
;
610 static void do_1pos(rvec xnew
, rvec xold
, rvec f
, real step
)
628 static void do_1pos3(rvec xnew
, rvec xold
, rvec f
, rvec step
)
646 static void directional_sd(rvec xold
[], rvec xnew
[], rvec acc_dir
[],
647 int start
, int homenr
, real step
)
651 for (i
= start
; i
< homenr
; i
++)
653 do_1pos(xnew
[i
], xold
[i
], acc_dir
[i
], step
);
657 static void shell_pos_sd(rvec xcur
[], rvec xnew
[], rvec f
[],
658 int ns
, t_shell s
[], int count
)
660 const real step_scale_min
= 0.8,
661 step_scale_increment
= 0.2,
662 step_scale_max
= 1.2,
663 step_scale_multiple
= (step_scale_max
- step_scale_min
) / step_scale_increment
;
667 real step_min
, step_max
;
672 for (i
= 0; (i
< ns
); i
++)
677 for (d
= 0; d
< DIM
; d
++)
679 s
[i
].step
[d
] = s
[i
].k_1
;
681 step_min
= min(step_min
, s
[i
].step
[d
]);
682 step_max
= max(step_max
, s
[i
].step
[d
]);
688 for (d
= 0; d
< DIM
; d
++)
690 dx
= xcur
[shell
][d
] - s
[i
].xold
[d
];
691 df
= f
[shell
][d
] - s
[i
].fold
[d
];
692 /* -dx/df gets used to generate an interpolated value, but would
693 * cause a NaN if df were binary-equal to zero. Values close to
694 * zero won't cause problems (because of the min() and max()), so
695 * just testing for binary inequality is OK. */
699 /* Scale the step size by a factor interpolated from
700 * step_scale_min to step_scale_max, as k_est goes from 0 to
701 * step_scale_multiple * s[i].step[d] */
703 step_scale_min
* s
[i
].step
[d
] +
704 step_scale_increment
* min(step_scale_multiple
* s
[i
].step
[d
], max(k_est
, 0));
709 if (gmx_numzero(dx
)) /* 0 == dx */
711 /* Likely this will never happen, but if it does just
712 * don't scale the step. */
716 s
[i
].step
[d
] *= step_scale_max
;
720 step_min
= min(step_min
, s
[i
].step
[d
]);
721 step_max
= max(step_max
, s
[i
].step
[d
]);
725 copy_rvec(xcur
[shell
], s
[i
].xold
);
726 copy_rvec(f
[shell
], s
[i
].fold
);
728 do_1pos3(xnew
[shell
], xcur
[shell
], f
[shell
], s
[i
].step
);
732 fprintf(debug
, "shell[%d] = %d\n", i
, shell
);
733 pr_rvec(debug
, 0, "fshell", f
[shell
], DIM
, TRUE
);
734 pr_rvec(debug
, 0, "xold", xcur
[shell
], DIM
, TRUE
);
735 pr_rvec(debug
, 0, "step", s
[i
].step
, DIM
, TRUE
);
736 pr_rvec(debug
, 0, "xnew", xnew
[shell
], DIM
, TRUE
);
740 printf("step %.3e %.3e\n", step_min
, step_max
);
744 static void decrease_step_size(int nshell
, t_shell s
[])
748 for (i
= 0; i
< nshell
; i
++)
750 svmul(0.8, s
[i
].step
, s
[i
].step
);
754 static void print_epot(FILE *fp
, gmx_int64_t mdstep
, int count
, real epot
, real df
,
755 int ndir
, real sf_dir
)
759 fprintf(fp
, "MDStep=%5s/%2d EPot: %12.8e, rmsF: %6.2e",
760 gmx_step_str(mdstep
, buf
), count
, epot
, df
);
763 fprintf(fp
, ", dir. rmsF: %6.2e\n", sqrt(sf_dir
/ndir
));
772 static real
rms_force(t_commrec
*cr
, rvec f
[], int ns
, t_shell s
[],
773 int ndir
, real
*sf_dir
, real
*Epot
)
779 for (i
= 0; i
< ns
; i
++)
782 buf
[0] += norm2(f
[shell
]);
791 gmx_sumd(4, buf
, cr
);
792 ntot
= (int)(buf
[1] + 0.5);
798 return (ntot
? sqrt(buf
[0]/ntot
) : 0);
801 static void check_pbc(FILE *fp
, rvec x
[], int shell
)
806 for (m
= 0; (m
< DIM
); m
++)
808 if (fabs(x
[shell
][m
]-x
[now
][m
]) > 0.3)
810 pr_rvecs(fp
, 0, "SHELL-X", x
+now
, 5);
816 static void dump_shells(FILE *fp
, rvec x
[], rvec f
[], real ftol
, int ns
, t_shell s
[])
823 for (i
= 0; (i
< ns
); i
++)
826 ff2
= iprod(f
[shell
], f
[shell
]);
829 fprintf(fp
, "SHELL %5d, force %10.5f %10.5f %10.5f, |f| %10.5f\n",
830 shell
, f
[shell
][XX
], f
[shell
][YY
], f
[shell
][ZZ
], sqrt(ff2
));
832 check_pbc(fp
, x
, shell
);
836 static void init_adir(FILE *log
, gmx_shellfc_t shfc
,
837 gmx_constr_t constr
, t_idef
*idef
, t_inputrec
*ir
,
838 t_commrec
*cr
, int dd_ac1
,
839 gmx_int64_t step
, t_mdatoms
*md
, int start
, int end
,
840 rvec
*x_old
, rvec
*x_init
, rvec
*x
,
841 rvec
*f
, rvec
*acc_dir
,
842 gmx_bool bMolPBC
, matrix box
,
843 real
*lambda
, real
*dvdlambda
, t_nrnb
*nrnb
)
850 unsigned short *ptype
;
853 if (DOMAINDECOMP(cr
))
861 if (n
> shfc
->adir_nalloc
)
863 shfc
->adir_nalloc
= over_alloc_dd(n
);
864 srenew(shfc
->adir_xnold
, shfc
->adir_nalloc
);
865 srenew(shfc
->adir_xnew
, shfc
->adir_nalloc
);
867 xnold
= shfc
->adir_xnold
;
868 xnew
= shfc
->adir_xnew
;
874 /* Does NOT work with freeze or acceleration groups (yet) */
875 for (n
= start
; n
< end
; n
++)
877 w_dt
= md
->invmass
[n
]*dt
;
879 for (d
= 0; d
< DIM
; d
++)
881 if ((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
))
883 xnold
[n
-start
][d
] = x
[n
][d
] - (x_init
[n
][d
] - x_old
[n
][d
]);
884 xnew
[n
-start
][d
] = 2*x
[n
][d
] - x_old
[n
][d
] + f
[n
][d
]*w_dt
*dt
;
888 xnold
[n
-start
][d
] = x
[n
][d
];
889 xnew
[n
-start
][d
] = x
[n
][d
];
893 constrain(log
, FALSE
, FALSE
, constr
, idef
, ir
, NULL
, cr
, step
, 0, md
,
894 x
, xnold
-start
, NULL
, bMolPBC
, box
,
895 lambda
[efptBONDED
], &(dvdlambda
[efptBONDED
]),
896 NULL
, NULL
, nrnb
, econqCoord
, FALSE
, 0, 0);
897 constrain(log
, FALSE
, FALSE
, constr
, idef
, ir
, NULL
, cr
, step
, 0, md
,
898 x
, xnew
-start
, NULL
, bMolPBC
, box
,
899 lambda
[efptBONDED
], &(dvdlambda
[efptBONDED
]),
900 NULL
, NULL
, nrnb
, econqCoord
, FALSE
, 0, 0);
902 for (n
= start
; n
< end
; n
++)
904 for (d
= 0; d
< DIM
; d
++)
907 -(2*x
[n
][d
]-xnold
[n
-start
][d
]-xnew
[n
-start
][d
])/sqr(dt
)
908 - f
[n
][d
]*md
->invmass
[n
];
910 clear_rvec(acc_dir
[n
]);
913 /* Project the acceleration on the old bond directions */
914 constrain(log
, FALSE
, FALSE
, constr
, idef
, ir
, NULL
, cr
, step
, 0, md
,
915 x_old
, xnew
-start
, acc_dir
, bMolPBC
, box
,
916 lambda
[efptBONDED
], &(dvdlambda
[efptBONDED
]),
917 NULL
, NULL
, nrnb
, econqDeriv_FlexCon
, FALSE
, 0, 0);
920 int relax_shell_flexcon(FILE *fplog
, t_commrec
*cr
, gmx_bool bVerbose
,
921 gmx_int64_t mdstep
, t_inputrec
*inputrec
,
922 gmx_bool bDoNS
, int force_flags
,
925 gmx_enerdata_t
*enerd
, t_fcdata
*fcd
,
926 t_state
*state
, rvec f
[],
929 t_nrnb
*nrnb
, gmx_wallcycle_t wcycle
,
931 gmx_groups_t
*groups
,
932 struct gmx_shellfc
*shfc
,
935 double t
, rvec mu_tot
,
936 gmx_bool
*bConverged
,
943 rvec
*pos
[2], *force
[2], *acc_dir
= NULL
, *x_old
= NULL
;
947 real ftol
, xiH
, xiS
, dum
= 0;
949 gmx_bool bCont
, bInit
;
950 int nat
, dd_ac0
, dd_ac1
= 0, i
;
951 int start
= 0, homenr
= md
->homenr
, end
= start
+homenr
, cg0
, cg1
;
952 int nflexcon
, g
, number_steps
, d
, Min
= 0, count
= 0;
953 #define Try (1-Min) /* At start Try = 1 */
955 bCont
= (mdstep
== inputrec
->init_step
) && inputrec
->bContinuation
;
956 bInit
= (mdstep
== inputrec
->init_step
) || shfc
->bRequireInit
;
957 ftol
= inputrec
->em_tol
;
958 number_steps
= inputrec
->niter
;
959 nshell
= shfc
->nshell
;
961 nflexcon
= shfc
->nflexcon
;
965 if (DOMAINDECOMP(cr
))
967 nat
= dd_natoms_vsite(cr
->dd
);
970 dd_get_constraint_range(cr
->dd
, &dd_ac0
, &dd_ac1
);
971 nat
= max(nat
, dd_ac1
);
979 if (nat
> shfc
->x_nalloc
)
981 /* Allocate local arrays */
982 shfc
->x_nalloc
= over_alloc_dd(nat
);
983 for (i
= 0; (i
< 2); i
++)
985 srenew(shfc
->x
[i
], shfc
->x_nalloc
);
986 srenew(shfc
->f
[i
], shfc
->x_nalloc
);
989 for (i
= 0; (i
< 2); i
++)
992 force
[i
] = shfc
->f
[i
];
995 if (bDoNS
&& inputrec
->ePBC
!= epbcNONE
&& !DOMAINDECOMP(cr
))
997 /* This is the only time where the coordinates are used
998 * before do_force is called, which normally puts all
999 * charge groups in the box.
1001 if (inputrec
->cutoff_scheme
== ecutsVERLET
)
1003 put_atoms_in_box_omp(fr
->ePBC
, state
->box
, md
->homenr
, state
->x
);
1009 put_charge_groups_in_box(fplog
, cg0
, cg1
, fr
->ePBC
, state
->box
,
1010 &(top
->cgs
), state
->x
, fr
->cg_cm
);
1015 mk_mshift(fplog
, graph
, fr
->ePBC
, state
->box
, state
->x
);
1019 /* After this all coordinate arrays will contain whole charge groups */
1022 shift_self(graph
, state
->box
, state
->x
);
1027 if (nat
> shfc
->flex_nalloc
)
1029 shfc
->flex_nalloc
= over_alloc_dd(nat
);
1030 srenew(shfc
->acc_dir
, shfc
->flex_nalloc
);
1031 srenew(shfc
->x_old
, shfc
->flex_nalloc
);
1033 acc_dir
= shfc
->acc_dir
;
1034 x_old
= shfc
->x_old
;
1035 for (i
= 0; i
< homenr
; i
++)
1037 for (d
= 0; d
< DIM
; d
++)
1040 state
->x
[start
+i
][d
] - state
->v
[start
+i
][d
]*inputrec
->delta_t
;
1045 /* Do a prediction of the shell positions */
1046 if (shfc
->bPredict
&& !bCont
)
1048 predict_shells(fplog
, state
->x
, state
->v
, inputrec
->delta_t
, nshell
, shell
,
1049 md
->massT
, NULL
, bInit
);
1052 /* do_force expected the charge groups to be in the box */
1055 unshift_self(graph
, state
->box
, state
->x
);
1058 /* Calculate the forces first time around */
1061 pr_rvecs(debug
, 0, "x b4 do_force", state
->x
+ start
, homenr
);
1063 do_force(fplog
, cr
, inputrec
, mdstep
, nrnb
, wcycle
, top
, groups
,
1064 state
->box
, state
->x
, &state
->hist
,
1065 force
[Min
], force_vir
, md
, enerd
, fcd
,
1066 state
->lambda
, graph
,
1067 fr
, vsite
, mu_tot
, t
, fp_field
, NULL
, bBornRadii
,
1068 (bDoNS
? GMX_FORCE_NS
: 0) | force_flags
);
1073 init_adir(fplog
, shfc
,
1074 constr
, idef
, inputrec
, cr
, dd_ac1
, mdstep
, md
, start
, end
,
1075 shfc
->x_old
-start
, state
->x
, state
->x
, force
[Min
],
1076 shfc
->acc_dir
-start
,
1077 fr
->bMolPBC
, state
->box
, state
->lambda
, &dum
, nrnb
);
1079 for (i
= start
; i
< end
; i
++)
1081 sf_dir
+= md
->massT
[i
]*norm2(shfc
->acc_dir
[i
-start
]);
1085 Epot
[Min
] = enerd
->term
[F_EPOT
];
1087 df
[Min
] = rms_force(cr
, shfc
->f
[Min
], nshell
, shell
, nflexcon
, &sf_dir
, &Epot
[Min
]);
1091 fprintf(debug
, "df = %g %g\n", df
[Min
], df
[Try
]);
1096 pr_rvecs(debug
, 0, "force0", force
[Min
], md
->nr
);
1099 if (nshell
+nflexcon
> 0)
1101 /* Copy x to pos[Min] & pos[Try]: during minimization only the
1102 * shell positions are updated, therefore the other particles must
1105 memcpy(pos
[Min
], state
->x
, nat
*sizeof(state
->x
[0]));
1106 memcpy(pos
[Try
], state
->x
, nat
*sizeof(state
->x
[0]));
1109 if (bVerbose
&& MASTER(cr
))
1111 print_epot(stdout
, mdstep
, 0, Epot
[Min
], df
[Min
], nflexcon
, sf_dir
);
1116 fprintf(debug
, "%17s: %14.10e\n",
1117 interaction_function
[F_EKIN
].longname
, enerd
->term
[F_EKIN
]);
1118 fprintf(debug
, "%17s: %14.10e\n",
1119 interaction_function
[F_EPOT
].longname
, enerd
->term
[F_EPOT
]);
1120 fprintf(debug
, "%17s: %14.10e\n",
1121 interaction_function
[F_ETOT
].longname
, enerd
->term
[F_ETOT
]);
1122 fprintf(debug
, "SHELLSTEP %s\n", gmx_step_str(mdstep
, sbuf
));
1125 /* First check whether we should do shells, or whether the force is
1126 * low enough even without minimization.
1128 *bConverged
= (df
[Min
] < ftol
);
1130 for (count
= 1; (!(*bConverged
) && (count
< number_steps
)); count
++)
1134 construct_vsites(vsite
, pos
[Min
], inputrec
->delta_t
, state
->v
,
1135 idef
->iparams
, idef
->il
,
1136 fr
->ePBC
, fr
->bMolPBC
, cr
, state
->box
);
1141 init_adir(fplog
, shfc
,
1142 constr
, idef
, inputrec
, cr
, dd_ac1
, mdstep
, md
, start
, end
,
1143 x_old
-start
, state
->x
, pos
[Min
], force
[Min
], acc_dir
-start
,
1144 fr
->bMolPBC
, state
->box
, state
->lambda
, &dum
, nrnb
);
1146 directional_sd(pos
[Min
], pos
[Try
], acc_dir
-start
, start
, end
,
1150 /* New positions, Steepest descent */
1151 shell_pos_sd(pos
[Min
], pos
[Try
], force
[Min
], nshell
, shell
, count
);
1153 /* do_force expected the charge groups to be in the box */
1156 unshift_self(graph
, state
->box
, pos
[Try
]);
1161 pr_rvecs(debug
, 0, "RELAX: pos[Min] ", pos
[Min
] + start
, homenr
);
1162 pr_rvecs(debug
, 0, "RELAX: pos[Try] ", pos
[Try
] + start
, homenr
);
1164 /* Try the new positions */
1165 do_force(fplog
, cr
, inputrec
, 1, nrnb
, wcycle
,
1166 top
, groups
, state
->box
, pos
[Try
], &state
->hist
,
1167 force
[Try
], force_vir
,
1168 md
, enerd
, fcd
, state
->lambda
, graph
,
1169 fr
, vsite
, mu_tot
, t
, fp_field
, NULL
, bBornRadii
,
1174 pr_rvecs(debug
, 0, "RELAX: force[Min]", force
[Min
] + start
, homenr
);
1175 pr_rvecs(debug
, 0, "RELAX: force[Try]", force
[Try
] + start
, homenr
);
1180 init_adir(fplog
, shfc
,
1181 constr
, idef
, inputrec
, cr
, dd_ac1
, mdstep
, md
, start
, end
,
1182 x_old
-start
, state
->x
, pos
[Try
], force
[Try
], acc_dir
-start
,
1183 fr
->bMolPBC
, state
->box
, state
->lambda
, &dum
, nrnb
);
1185 for (i
= start
; i
< end
; i
++)
1187 sf_dir
+= md
->massT
[i
]*norm2(acc_dir
[i
-start
]);
1191 Epot
[Try
] = enerd
->term
[F_EPOT
];
1193 df
[Try
] = rms_force(cr
, force
[Try
], nshell
, shell
, nflexcon
, &sf_dir
, &Epot
[Try
]);
1197 fprintf(debug
, "df = %g %g\n", df
[Min
], df
[Try
]);
1204 pr_rvecs(debug
, 0, "F na do_force", force
[Try
] + start
, homenr
);
1208 fprintf(debug
, "SHELL ITER %d\n", count
);
1209 dump_shells(debug
, pos
[Try
], force
[Try
], ftol
, nshell
, shell
);
1213 if (bVerbose
&& MASTER(cr
))
1215 print_epot(stdout
, mdstep
, count
, Epot
[Try
], df
[Try
], nflexcon
, sf_dir
);
1218 *bConverged
= (df
[Try
] < ftol
);
1220 if ((df
[Try
] < df
[Min
]))
1224 fprintf(debug
, "Swapping Min and Try\n");
1228 /* Correct the velocities for the flexible constraints */
1229 invdt
= 1/inputrec
->delta_t
;
1230 for (i
= start
; i
< end
; i
++)
1232 for (d
= 0; d
< DIM
; d
++)
1234 state
->v
[i
][d
] += (pos
[Try
][i
][d
] - pos
[Min
][i
][d
])*invdt
;
1242 decrease_step_size(nshell
, shell
);
1245 if (MASTER(cr
) && !(*bConverged
))
1247 /* Note that the energies and virial are incorrect when not converged */
1251 "step %s: EM did not converge in %d iterations, RMS force %.3f\n",
1252 gmx_step_str(mdstep
, sbuf
), number_steps
, df
[Min
]);
1255 "step %s: EM did not converge in %d iterations, RMS force %.3f\n",
1256 gmx_step_str(mdstep
, sbuf
), number_steps
, df
[Min
]);
1259 /* Copy back the coordinates and the forces */
1260 memcpy(state
->x
, pos
[Min
], nat
*sizeof(state
->x
[0]));
1261 memcpy(f
, force
[Min
], nat
*sizeof(f
[0]));