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-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015, 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.
39 * \brief This file defines integrators for energy minimization
41 * \author Berk Hess <hess@kth.se>
42 * \author Erik Lindahl <erik@kth.se>
43 * \ingroup module_mdlib
57 #include "gromacs/commandline/filenm.h"
58 #include "gromacs/domdec/domdec.h"
59 #include "gromacs/domdec/domdec_struct.h"
60 #include "gromacs/ewald/pme.h"
61 #include "gromacs/fileio/confio.h"
62 #include "gromacs/fileio/mtxio.h"
63 #include "gromacs/gmxlib/gmx_omp_nthreads.h"
64 #include "gromacs/gmxlib/md_logging.h"
65 #include "gromacs/gmxlib/network.h"
66 #include "gromacs/gmxlib/nrnb.h"
67 #include "gromacs/imd/imd.h"
68 #include "gromacs/linearalgebra/sparsematrix.h"
69 #include "gromacs/listed-forces/manage-threading.h"
70 #include "gromacs/math/vec.h"
71 #include "gromacs/mdlib/constr.h"
72 #include "gromacs/mdlib/force.h"
73 #include "gromacs/mdlib/forcerec.h"
74 #include "gromacs/mdlib/md_support.h"
75 #include "gromacs/mdlib/mdatoms.h"
76 #include "gromacs/mdlib/mdebin.h"
77 #include "gromacs/mdlib/mdrun.h"
78 #include "gromacs/mdlib/ns.h"
79 #include "gromacs/mdlib/sim_util.h"
80 #include "gromacs/mdlib/tgroup.h"
81 #include "gromacs/mdlib/trajectory_writing.h"
82 #include "gromacs/mdlib/update.h"
83 #include "gromacs/mdlib/vsite.h"
84 #include "gromacs/mdtypes/commrec.h"
85 #include "gromacs/mdtypes/inputrec.h"
86 #include "gromacs/mdtypes/md_enums.h"
87 #include "gromacs/pbcutil/mshift.h"
88 #include "gromacs/pbcutil/pbc.h"
89 #include "gromacs/timing/wallcycle.h"
90 #include "gromacs/timing/walltime_accounting.h"
91 #include "gromacs/topology/mtop_util.h"
92 #include "gromacs/utility/cstringutil.h"
93 #include "gromacs/utility/exceptions.h"
94 #include "gromacs/utility/fatalerror.h"
95 #include "gromacs/utility/smalloc.h"
97 //! Utility structure for manipulating states during EM
99 //! Copy of the global state
105 //! Norm of the force
113 //! Initiate em_state_t structure and return pointer to it
114 static em_state_t
*init_em_state()
120 /* does this need to be here? Should the array be declared differently (staticaly)in the state definition? */
121 snew(ems
->s
.lambda
, efptNR
);
126 //! Print the EM starting conditions
127 static void print_em_start(FILE *fplog
,
129 gmx_walltime_accounting_t walltime_accounting
,
130 gmx_wallcycle_t wcycle
,
133 walltime_accounting_start(walltime_accounting
);
134 wallcycle_start(wcycle
, ewcRUN
);
135 print_start(fplog
, cr
, walltime_accounting
, name
);
138 //! Stop counting time for EM
139 static void em_time_end(gmx_walltime_accounting_t walltime_accounting
,
140 gmx_wallcycle_t wcycle
)
142 wallcycle_stop(wcycle
, ewcRUN
);
144 walltime_accounting_end(walltime_accounting
);
147 //! Printing a log file and console header
148 static void sp_header(FILE *out
, const char *minimizer
, real ftol
, int nsteps
)
151 fprintf(out
, "%s:\n", minimizer
);
152 fprintf(out
, " Tolerance (Fmax) = %12.5e\n", ftol
);
153 fprintf(out
, " Number of steps = %12d\n", nsteps
);
156 //! Print warning message
157 static void warn_step(FILE *fp
, real ftol
, gmx_bool bLastStep
, gmx_bool bConstrain
)
163 "\nEnergy minimization reached the maximum number "
164 "of steps before the forces reached the requested "
165 "precision Fmax < %g.\n", ftol
);
170 "\nEnergy minimization has stopped, but the forces have "
171 "not converged to the requested precision Fmax < %g (which "
172 "may not be possible for your system). It stopped "
173 "because the algorithm tried to make a new step whose size "
174 "was too small, or there was no change in the energy since "
175 "last step. Either way, we regard the minimization as "
176 "converged to within the available machine precision, "
177 "given your starting configuration and EM parameters.\n%s%s",
179 sizeof(real
) < sizeof(double) ?
180 "\nDouble precision normally gives you higher accuracy, but "
181 "this is often not needed for preparing to run molecular "
185 "You might need to increase your constraint accuracy, or turn\n"
186 "off constraints altogether (set constraints = none in mdp file)\n" :
189 fputs(wrap_lines(buffer
, 78, 0, FALSE
), fp
);
192 //! Print message about convergence of the EM
193 static void print_converged(FILE *fp
, const char *alg
, real ftol
,
194 gmx_int64_t count
, gmx_bool bDone
, gmx_int64_t nsteps
,
195 real epot
, real fmax
, int nfmax
, real fnorm
)
197 char buf
[STEPSTRSIZE
];
201 fprintf(fp
, "\n%s converged to Fmax < %g in %s steps\n",
202 alg
, ftol
, gmx_step_str(count
, buf
));
204 else if (count
< nsteps
)
206 fprintf(fp
, "\n%s converged to machine precision in %s steps,\n"
207 "but did not reach the requested Fmax < %g.\n",
208 alg
, gmx_step_str(count
, buf
), ftol
);
212 fprintf(fp
, "\n%s did not converge to Fmax < %g in %s steps.\n",
213 alg
, ftol
, gmx_step_str(count
, buf
));
217 fprintf(fp
, "Potential Energy = %21.14e\n", epot
);
218 fprintf(fp
, "Maximum force = %21.14e on atom %d\n", fmax
, nfmax
+1);
219 fprintf(fp
, "Norm of force = %21.14e\n", fnorm
);
221 fprintf(fp
, "Potential Energy = %14.7e\n", epot
);
222 fprintf(fp
, "Maximum force = %14.7e on atom %d\n", fmax
, nfmax
+1);
223 fprintf(fp
, "Norm of force = %14.7e\n", fnorm
);
227 //! Compute the norm and max of the force array in parallel
228 static void get_f_norm_max(t_commrec
*cr
,
229 t_grpopts
*opts
, t_mdatoms
*mdatoms
, rvec
*f
,
230 real
*fnorm
, real
*fmax
, int *a_fmax
)
234 int la_max
, a_max
, start
, end
, i
, m
, gf
;
236 /* This routine finds the largest force and returns it.
237 * On parallel machines the global max is taken.
243 end
= mdatoms
->homenr
;
244 if (mdatoms
->cFREEZE
)
246 for (i
= start
; i
< end
; i
++)
248 gf
= mdatoms
->cFREEZE
[i
];
250 for (m
= 0; m
< DIM
; m
++)
252 if (!opts
->nFreeze
[gf
][m
])
267 for (i
= start
; i
< end
; i
++)
279 if (la_max
>= 0 && DOMAINDECOMP(cr
))
281 a_max
= cr
->dd
->gatindex
[la_max
];
289 snew(sum
, 2*cr
->nnodes
+1);
290 sum
[2*cr
->nodeid
] = fmax2
;
291 sum
[2*cr
->nodeid
+1] = a_max
;
292 sum
[2*cr
->nnodes
] = fnorm2
;
293 gmx_sumd(2*cr
->nnodes
+1, sum
, cr
);
294 fnorm2
= sum
[2*cr
->nnodes
];
295 /* Determine the global maximum */
296 for (i
= 0; i
< cr
->nnodes
; i
++)
298 if (sum
[2*i
] > fmax2
)
301 a_max
= (int)(sum
[2*i
+1] + 0.5);
309 *fnorm
= sqrt(fnorm2
);
321 //! Compute the norm of the force
322 static void get_state_f_norm_max(t_commrec
*cr
,
323 t_grpopts
*opts
, t_mdatoms
*mdatoms
,
326 get_f_norm_max(cr
, opts
, mdatoms
, ems
->f
, &ems
->fnorm
, &ems
->fmax
, &ems
->a_fmax
);
329 //! Initialize the energy minimization
330 void init_em(FILE *fplog
, const char *title
,
331 t_commrec
*cr
, t_inputrec
*ir
,
332 t_state
*state_global
, gmx_mtop_t
*top_global
,
333 em_state_t
*ems
, gmx_localtop_t
**top
,
334 rvec
**f
, rvec
**f_global
,
335 t_nrnb
*nrnb
, rvec mu_tot
,
336 t_forcerec
*fr
, gmx_enerdata_t
**enerd
,
337 t_graph
**graph
, t_mdatoms
*mdatoms
, gmx_global_stat_t
*gstat
,
338 gmx_vsite_t
*vsite
, gmx_constr_t constr
,
339 int nfile
, const t_filenm fnm
[],
340 gmx_mdoutf_t
*outf
, t_mdebin
**mdebin
,
341 int imdport
, unsigned long gmx_unused Flags
,
342 gmx_wallcycle_t wcycle
)
349 fprintf(fplog
, "Initiating %s\n", title
);
352 state_global
->ngtc
= 0;
354 /* Initialize lambda variables */
355 initialize_lambdas(fplog
, ir
, &(state_global
->fep_state
), state_global
->lambda
, NULL
);
359 /* Interactive molecular dynamics */
360 init_IMD(ir
, cr
, top_global
, fplog
, 1, state_global
->x
,
361 nfile
, fnm
, NULL
, imdport
, Flags
);
363 if (DOMAINDECOMP(cr
))
365 *top
= dd_init_local_top(top_global
);
367 dd_init_local_state(cr
->dd
, state_global
, &ems
->s
);
371 /* Distribute the charge groups over the nodes from the master node */
372 dd_partition_system(fplog
, ir
->init_step
, cr
, TRUE
, 1,
373 state_global
, top_global
, ir
,
374 &ems
->s
, &ems
->f
, mdatoms
, *top
,
375 fr
, vsite
, NULL
, constr
,
377 dd_store_state(cr
->dd
, &ems
->s
);
381 snew(*f_global
, top_global
->natoms
);
391 snew(*f
, top_global
->natoms
);
393 /* Just copy the state */
394 ems
->s
= *state_global
;
395 snew(ems
->s
.x
, ems
->s
.nalloc
);
396 snew(ems
->f
, ems
->s
.nalloc
);
397 for (i
= 0; i
< state_global
->natoms
; i
++)
399 copy_rvec(state_global
->x
[i
], ems
->s
.x
[i
]);
401 copy_mat(state_global
->box
, ems
->s
.box
);
403 *top
= gmx_mtop_generate_local_top(top_global
, ir
->efep
!= efepNO
);
406 forcerec_set_excl_load(fr
, *top
);
408 setup_bonded_threading(fr
, &(*top
)->idef
);
410 if (ir
->ePBC
!= epbcNONE
&& !fr
->bMolPBC
)
412 *graph
= mk_graph(fplog
, &((*top
)->idef
), 0, top_global
->natoms
, FALSE
, FALSE
);
419 atoms2md(top_global
, ir
, 0, NULL
, top_global
->natoms
, mdatoms
);
420 update_mdatoms(mdatoms
, state_global
->lambda
[efptFEP
]);
424 set_vsite_top(vsite
, *top
, mdatoms
, cr
);
430 if (ir
->eConstrAlg
== econtSHAKE
&&
431 gmx_mtop_ftype_count(top_global
, F_CONSTR
) > 0)
433 gmx_fatal(FARGS
, "Can not do energy minimization with %s, use %s\n",
434 econstr_names
[econtSHAKE
], econstr_names
[econtLINCS
]);
437 if (!DOMAINDECOMP(cr
))
439 set_constraints(constr
, *top
, ir
, mdatoms
, cr
);
442 if (!ir
->bContinuation
)
444 /* Constrain the starting coordinates */
446 constrain(PAR(cr
) ? NULL
: fplog
, TRUE
, TRUE
, constr
, &(*top
)->idef
,
447 ir
, cr
, -1, 0, 1.0, mdatoms
,
448 ems
->s
.x
, ems
->s
.x
, NULL
, fr
->bMolPBC
, ems
->s
.box
,
449 ems
->s
.lambda
[efptFEP
], &dvdl_constr
,
450 NULL
, NULL
, nrnb
, econqCoord
);
456 *gstat
= global_stat_init(ir
);
463 *outf
= init_mdoutf(fplog
, nfile
, fnm
, 0, cr
, ir
, top_global
, NULL
, wcycle
);
466 init_enerdata(top_global
->groups
.grps
[egcENER
].nr
, ir
->fepvals
->n_lambda
,
471 /* Init bin for energy stuff */
472 *mdebin
= init_mdebin(mdoutf_get_fp_ene(*outf
), top_global
, ir
, NULL
);
476 calc_shifts(ems
->s
.box
, fr
->shift_vec
);
479 //! Finalize the minimization
480 static void finish_em(t_commrec
*cr
, gmx_mdoutf_t outf
,
481 gmx_walltime_accounting_t walltime_accounting
,
482 gmx_wallcycle_t wcycle
)
484 if (!(cr
->duty
& DUTY_PME
))
486 /* Tell the PME only node to finish */
487 gmx_pme_send_finish(cr
);
492 em_time_end(walltime_accounting
, wcycle
);
495 //! Swap two different EM states during minimization
496 static void swap_em_state(em_state_t
*ems1
, em_state_t
*ems2
)
505 //! Copy coordinate from an EM state to a "normal" state structure
506 static void copy_em_coords(em_state_t
*ems
, t_state
*state
)
510 for (i
= 0; (i
< state
->natoms
); i
++)
512 copy_rvec(ems
->s
.x
[i
], state
->x
[i
]);
516 //! Save the EM trajectory
517 static void write_em_traj(FILE *fplog
, t_commrec
*cr
,
519 gmx_bool bX
, gmx_bool bF
, const char *confout
,
520 gmx_mtop_t
*top_global
,
521 t_inputrec
*ir
, gmx_int64_t step
,
523 t_state
*state_global
, rvec
*f_global
)
526 gmx_bool bIMDout
= FALSE
;
529 /* Shall we do IMD output? */
532 bIMDout
= do_per_step(step
, IMD_get_step(ir
->imd
->setup
));
535 if ((bX
|| bF
|| bIMDout
|| confout
!= NULL
) && !DOMAINDECOMP(cr
))
537 copy_em_coords(state
, state_global
);
544 mdof_flags
|= MDOF_X
;
548 mdof_flags
|= MDOF_F
;
551 /* If we want IMD output, set appropriate MDOF flag */
554 mdof_flags
|= MDOF_IMD
;
557 mdoutf_write_to_trajectory_files(fplog
, cr
, outf
, mdof_flags
,
558 top_global
, step
, (double)step
,
559 &state
->s
, state_global
, state
->f
, f_global
);
561 if (confout
!= NULL
&& MASTER(cr
))
563 if (ir
->ePBC
!= epbcNONE
&& !ir
->bPeriodicMols
&& DOMAINDECOMP(cr
))
565 /* Make molecules whole only for confout writing */
566 do_pbc_mtop(fplog
, ir
->ePBC
, state_global
->box
, top_global
,
570 write_sto_conf_mtop(confout
,
571 *top_global
->name
, top_global
,
572 state_global
->x
, NULL
, ir
->ePBC
, state_global
->box
);
576 //! Do one minimization step
577 static void do_em_step(t_commrec
*cr
, t_inputrec
*ir
, t_mdatoms
*md
,
579 em_state_t
*ems1
, real a
, rvec
*f
, em_state_t
*ems2
,
580 gmx_constr_t constr
, gmx_localtop_t
*top
,
581 t_nrnb
*nrnb
, gmx_wallcycle_t wcycle
,
590 int nthreads gmx_unused
;
595 if (DOMAINDECOMP(cr
) && s1
->ddp_count
!= cr
->dd
->ddp_count
)
597 gmx_incons("state mismatch in do_em_step");
600 s2
->flags
= s1
->flags
;
602 if (s2
->nalloc
!= s1
->nalloc
)
604 s2
->nalloc
= s1
->nalloc
;
605 srenew(s2
->x
, s1
->nalloc
);
606 srenew(ems2
->f
, s1
->nalloc
);
607 if (s2
->flags
& (1<<estCGP
))
609 srenew(s2
->cg_p
, s1
->nalloc
);
613 s2
->natoms
= s1
->natoms
;
614 copy_mat(s1
->box
, s2
->box
);
615 /* Copy free energy state */
616 for (i
= 0; i
< efptNR
; i
++)
618 s2
->lambda
[i
] = s1
->lambda
[i
];
620 copy_mat(s1
->box
, s2
->box
);
628 // cppcheck-suppress unreadVariable
629 nthreads
= gmx_omp_nthreads_get(emntUpdate
);
630 #pragma omp parallel num_threads(nthreads)
635 #pragma omp for schedule(static) nowait
636 for (i
= start
; i
< end
; i
++)
644 for (m
= 0; m
< DIM
; m
++)
646 if (ir
->opts
.nFreeze
[gf
][m
])
652 x2
[i
][m
] = x1
[i
][m
] + a
*f
[i
][m
];
656 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
659 if (s2
->flags
& (1<<estCGP
))
661 /* Copy the CG p vector */
664 #pragma omp for schedule(static) nowait
665 for (i
= start
; i
< end
; i
++)
667 // Trivial OpenMP block that does not throw
668 copy_rvec(x1
[i
], x2
[i
]);
672 if (DOMAINDECOMP(cr
))
674 s2
->ddp_count
= s1
->ddp_count
;
675 if (s2
->cg_gl_nalloc
< s1
->cg_gl_nalloc
)
678 s2
->cg_gl_nalloc
= s1
->cg_gl_nalloc
;
681 srenew(s2
->cg_gl
, s2
->cg_gl_nalloc
);
683 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
686 s2
->ncg_gl
= s1
->ncg_gl
;
687 #pragma omp for schedule(static) nowait
688 for (i
= 0; i
< s2
->ncg_gl
; i
++)
690 s2
->cg_gl
[i
] = s1
->cg_gl
[i
];
692 s2
->ddp_count_cg_gl
= s1
->ddp_count_cg_gl
;
698 wallcycle_start(wcycle
, ewcCONSTR
);
700 constrain(NULL
, TRUE
, TRUE
, constr
, &top
->idef
,
701 ir
, cr
, count
, 0, 1.0, md
,
702 s1
->x
, s2
->x
, NULL
, bMolPBC
, s2
->box
,
703 s2
->lambda
[efptBONDED
], &dvdl_constr
,
704 NULL
, NULL
, nrnb
, econqCoord
);
705 wallcycle_stop(wcycle
, ewcCONSTR
);
709 //! Prepare EM for using domain decomposition parallellization
710 static void em_dd_partition_system(FILE *fplog
, int step
, t_commrec
*cr
,
711 gmx_mtop_t
*top_global
, t_inputrec
*ir
,
712 em_state_t
*ems
, gmx_localtop_t
*top
,
713 t_mdatoms
*mdatoms
, t_forcerec
*fr
,
714 gmx_vsite_t
*vsite
, gmx_constr_t constr
,
715 t_nrnb
*nrnb
, gmx_wallcycle_t wcycle
)
717 /* Repartition the domain decomposition */
718 dd_partition_system(fplog
, step
, cr
, FALSE
, 1,
719 NULL
, top_global
, ir
,
721 mdatoms
, top
, fr
, vsite
, NULL
, constr
,
722 nrnb
, wcycle
, FALSE
);
723 dd_store_state(cr
->dd
, &ems
->s
);
726 //! De one energy evaluation
727 static void evaluate_energy(FILE *fplog
, t_commrec
*cr
,
728 gmx_mtop_t
*top_global
,
729 em_state_t
*ems
, gmx_localtop_t
*top
,
730 t_inputrec
*inputrec
,
731 t_nrnb
*nrnb
, gmx_wallcycle_t wcycle
,
732 gmx_global_stat_t gstat
,
733 gmx_vsite_t
*vsite
, gmx_constr_t constr
,
735 t_graph
*graph
, t_mdatoms
*mdatoms
,
736 t_forcerec
*fr
, rvec mu_tot
,
737 gmx_enerdata_t
*enerd
, tensor vir
, tensor pres
,
738 gmx_int64_t count
, gmx_bool bFirst
)
742 tensor force_vir
, shake_vir
, ekin
;
743 real dvdl_constr
, prescorr
, enercorr
, dvdlcorr
;
746 /* Set the time to the initial time, the time does not change during EM */
747 t
= inputrec
->init_t
;
750 (DOMAINDECOMP(cr
) && ems
->s
.ddp_count
< cr
->dd
->ddp_count
))
752 /* This is the first state or an old state used before the last ns */
758 if (inputrec
->nstlist
> 0)
766 construct_vsites(vsite
, ems
->s
.x
, 1, NULL
,
767 top
->idef
.iparams
, top
->idef
.il
,
768 fr
->ePBC
, fr
->bMolPBC
, cr
, ems
->s
.box
);
771 if (DOMAINDECOMP(cr
) && bNS
)
773 /* Repartition the domain decomposition */
774 em_dd_partition_system(fplog
, count
, cr
, top_global
, inputrec
,
775 ems
, top
, mdatoms
, fr
, vsite
, constr
,
779 /* Calc force & energy on new trial position */
780 /* do_force always puts the charge groups in the box and shifts again
781 * We do not unshift, so molecules are always whole in congrad.c
783 do_force(fplog
, cr
, inputrec
,
784 count
, nrnb
, wcycle
, top
, &top_global
->groups
,
785 ems
->s
.box
, ems
->s
.x
, &ems
->s
.hist
,
786 ems
->f
, force_vir
, mdatoms
, enerd
, fcd
,
787 ems
->s
.lambda
, graph
, fr
, vsite
, mu_tot
, t
, NULL
, NULL
, TRUE
,
788 GMX_FORCE_STATECHANGED
| GMX_FORCE_ALLFORCES
|
789 GMX_FORCE_VIRIAL
| GMX_FORCE_ENERGY
|
790 (bNS
? GMX_FORCE_NS
: 0));
792 /* Clear the unused shake virial and pressure */
793 clear_mat(shake_vir
);
796 /* Communicate stuff when parallel */
797 if (PAR(cr
) && inputrec
->eI
!= eiNM
)
799 wallcycle_start(wcycle
, ewcMoveE
);
801 global_stat(fplog
, gstat
, cr
, enerd
, force_vir
, shake_vir
, mu_tot
,
802 inputrec
, NULL
, NULL
, NULL
, 1, &terminate
,
803 top_global
, &ems
->s
, FALSE
,
808 wallcycle_stop(wcycle
, ewcMoveE
);
811 /* Calculate long range corrections to pressure and energy */
812 calc_dispcorr(inputrec
, fr
, top_global
->natoms
, ems
->s
.box
, ems
->s
.lambda
[efptVDW
],
813 pres
, force_vir
, &prescorr
, &enercorr
, &dvdlcorr
);
814 enerd
->term
[F_DISPCORR
] = enercorr
;
815 enerd
->term
[F_EPOT
] += enercorr
;
816 enerd
->term
[F_PRES
] += prescorr
;
817 enerd
->term
[F_DVDL
] += dvdlcorr
;
819 ems
->epot
= enerd
->term
[F_EPOT
];
823 /* Project out the constraint components of the force */
824 wallcycle_start(wcycle
, ewcCONSTR
);
826 constrain(NULL
, FALSE
, FALSE
, constr
, &top
->idef
,
827 inputrec
, cr
, count
, 0, 1.0, mdatoms
,
828 ems
->s
.x
, ems
->f
, ems
->f
, fr
->bMolPBC
, ems
->s
.box
,
829 ems
->s
.lambda
[efptBONDED
], &dvdl_constr
,
830 NULL
, &shake_vir
, nrnb
, econqForceDispl
);
831 enerd
->term
[F_DVDL_CONSTR
] += dvdl_constr
;
832 m_add(force_vir
, shake_vir
, vir
);
833 wallcycle_stop(wcycle
, ewcCONSTR
);
837 copy_mat(force_vir
, vir
);
841 enerd
->term
[F_PRES
] =
842 calc_pres(fr
->ePBC
, inputrec
->nwall
, ems
->s
.box
, ekin
, vir
, pres
);
844 sum_dhdl(enerd
, ems
->s
.lambda
, inputrec
->fepvals
);
846 if (EI_ENERGY_MINIMIZATION(inputrec
->eI
))
848 get_state_f_norm_max(cr
, &(inputrec
->opts
), mdatoms
, ems
);
852 //! Parallel utility summing energies and forces
853 static double reorder_partsum(t_commrec
*cr
, t_grpopts
*opts
, t_mdatoms
*mdatoms
,
854 gmx_mtop_t
*top_global
,
855 em_state_t
*s_min
, em_state_t
*s_b
)
859 int ncg
, *cg_gl
, *index
, c
, cg
, i
, a0
, a1
, a
, gf
, m
;
861 unsigned char *grpnrFREEZE
;
865 fprintf(debug
, "Doing reorder_partsum\n");
871 cgs_gl
= dd_charge_groups_global(cr
->dd
);
872 index
= cgs_gl
->index
;
874 /* Collect fm in a global vector fmg.
875 * This conflicts with the spirit of domain decomposition,
876 * but to fully optimize this a much more complicated algorithm is required.
878 snew(fmg
, top_global
->natoms
);
880 ncg
= s_min
->s
.ncg_gl
;
881 cg_gl
= s_min
->s
.cg_gl
;
883 for (c
= 0; c
< ncg
; c
++)
888 for (a
= a0
; a
< a1
; a
++)
890 copy_rvec(fm
[i
], fmg
[a
]);
894 gmx_sum(top_global
->natoms
*3, fmg
[0], cr
);
896 /* Now we will determine the part of the sum for the cgs in state s_b */
898 cg_gl
= s_b
->s
.cg_gl
;
902 grpnrFREEZE
= top_global
->groups
.grpnr
[egcFREEZE
];
903 for (c
= 0; c
< ncg
; c
++)
908 for (a
= a0
; a
< a1
; a
++)
910 if (mdatoms
->cFREEZE
&& grpnrFREEZE
)
914 for (m
= 0; m
< DIM
; m
++)
916 if (!opts
->nFreeze
[gf
][m
])
918 partsum
+= (fb
[i
][m
] - fmg
[a
][m
])*fb
[i
][m
];
930 //! Print some stuff, like beta, whatever that means.
931 static real
pr_beta(t_commrec
*cr
, t_grpopts
*opts
, t_mdatoms
*mdatoms
,
932 gmx_mtop_t
*top_global
,
933 em_state_t
*s_min
, em_state_t
*s_b
)
939 /* This is just the classical Polak-Ribiere calculation of beta;
940 * it looks a bit complicated since we take freeze groups into account,
941 * and might have to sum it in parallel runs.
944 if (!DOMAINDECOMP(cr
) ||
945 (s_min
->s
.ddp_count
== cr
->dd
->ddp_count
&&
946 s_b
->s
.ddp_count
== cr
->dd
->ddp_count
))
952 /* This part of code can be incorrect with DD,
953 * since the atom ordering in s_b and s_min might differ.
955 for (i
= 0; i
< mdatoms
->homenr
; i
++)
957 if (mdatoms
->cFREEZE
)
959 gf
= mdatoms
->cFREEZE
[i
];
961 for (m
= 0; m
< DIM
; m
++)
963 if (!opts
->nFreeze
[gf
][m
])
965 sum
+= (fb
[i
][m
] - fm
[i
][m
])*fb
[i
][m
];
972 /* We need to reorder cgs while summing */
973 sum
= reorder_partsum(cr
, opts
, mdatoms
, top_global
, s_min
, s_b
);
977 gmx_sumd(1, &sum
, cr
);
980 return sum
/sqr(s_min
->fnorm
);
986 /*! \brief Do conjugate gradients minimization
987 \copydoc integrator_t(FILE *fplog, t_commrec *cr,
988 int nfile, const t_filenm fnm[],
989 const gmx_output_env_t *oenv, gmx_bool bVerbose,
991 gmx_vsite_t *vsite, gmx_constr_t constr,
993 t_inputrec *inputrec,
994 gmx_mtop_t *top_global, t_fcdata *fcd,
995 t_state *state_global,
997 t_nrnb *nrnb, gmx_wallcycle_t wcycle,
1000 int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
1001 gmx_membed_t *membed,
1002 real cpt_period, real max_hours,
1004 unsigned long Flags,
1005 gmx_walltime_accounting_t walltime_accounting)
1007 double do_cg(FILE *fplog
, t_commrec
*cr
,
1008 int nfile
, const t_filenm fnm
[],
1009 const gmx_output_env_t gmx_unused
*oenv
, gmx_bool bVerbose
,
1010 int gmx_unused nstglobalcomm
,
1011 gmx_vsite_t
*vsite
, gmx_constr_t constr
,
1012 int gmx_unused stepout
,
1013 t_inputrec
*inputrec
,
1014 gmx_mtop_t
*top_global
, t_fcdata
*fcd
,
1015 t_state
*state_global
,
1017 t_nrnb
*nrnb
, gmx_wallcycle_t wcycle
,
1018 gmx_edsam_t gmx_unused ed
,
1020 int gmx_unused repl_ex_nst
, int gmx_unused repl_ex_nex
, int gmx_unused repl_ex_seed
,
1021 gmx_membed_t
* /*membed*/,
1022 real gmx_unused cpt_period
, real gmx_unused max_hours
,
1024 unsigned long gmx_unused Flags
,
1025 gmx_walltime_accounting_t walltime_accounting
)
1027 const char *CG
= "Polak-Ribiere Conjugate Gradients";
1029 em_state_t
*s_min
, *s_a
, *s_b
, *s_c
;
1030 gmx_localtop_t
*top
;
1031 gmx_enerdata_t
*enerd
;
1033 gmx_global_stat_t gstat
;
1035 rvec
*f_global
, *p
, *sf
;
1036 double gpa
, gpb
, gpc
, tmp
, minstep
;
1039 real a
, b
, c
, beta
= 0.0;
1043 gmx_bool converged
, foundlower
;
1045 gmx_bool do_log
= FALSE
, do_ene
= FALSE
, do_x
, do_f
;
1047 int number_steps
, neval
= 0, nstcg
= inputrec
->nstcgsteep
;
1049 int i
, m
, gf
, step
, nminstep
;
1053 s_min
= init_em_state();
1054 s_a
= init_em_state();
1055 s_b
= init_em_state();
1056 s_c
= init_em_state();
1058 /* Init em and store the local state in s_min */
1059 init_em(fplog
, CG
, cr
, inputrec
,
1060 state_global
, top_global
, s_min
, &top
, &f
, &f_global
,
1061 nrnb
, mu_tot
, fr
, &enerd
, &graph
, mdatoms
, &gstat
, vsite
, constr
,
1062 nfile
, fnm
, &outf
, &mdebin
, imdport
, Flags
, wcycle
);
1064 /* Print to log file */
1065 print_em_start(fplog
, cr
, walltime_accounting
, wcycle
, CG
);
1067 /* Max number of steps */
1068 number_steps
= inputrec
->nsteps
;
1072 sp_header(stderr
, CG
, inputrec
->em_tol
, number_steps
);
1076 sp_header(fplog
, CG
, inputrec
->em_tol
, number_steps
);
1079 /* Call the force routine and some auxiliary (neighboursearching etc.) */
1080 /* do_force always puts the charge groups in the box and shifts again
1081 * We do not unshift, so molecules are always whole in congrad.c
1083 evaluate_energy(fplog
, cr
,
1084 top_global
, s_min
, top
,
1085 inputrec
, nrnb
, wcycle
, gstat
,
1086 vsite
, constr
, fcd
, graph
, mdatoms
, fr
,
1087 mu_tot
, enerd
, vir
, pres
, -1, TRUE
);
1092 /* Copy stuff to the energy bin for easy printing etc. */
1093 upd_mdebin(mdebin
, FALSE
, FALSE
, (double)step
,
1094 mdatoms
->tmass
, enerd
, &s_min
->s
, inputrec
->fepvals
, inputrec
->expandedvals
, s_min
->s
.box
,
1095 NULL
, NULL
, vir
, pres
, NULL
, mu_tot
, constr
);
1097 print_ebin_header(fplog
, step
, step
);
1098 print_ebin(mdoutf_get_fp_ene(outf
), TRUE
, FALSE
, FALSE
, fplog
, step
, step
, eprNORMAL
,
1099 mdebin
, fcd
, &(top_global
->groups
), &(inputrec
->opts
));
1103 /* Estimate/guess the initial stepsize */
1104 stepsize
= inputrec
->em_stepsize
/s_min
->fnorm
;
1108 double sqrtNumAtoms
= sqrt(static_cast<double>(state_global
->natoms
));
1109 fprintf(stderr
, " F-max = %12.5e on atom %d\n",
1110 s_min
->fmax
, s_min
->a_fmax
+1);
1111 fprintf(stderr
, " F-Norm = %12.5e\n",
1112 s_min
->fnorm
/sqrtNumAtoms
);
1113 fprintf(stderr
, "\n");
1114 /* and copy to the log file too... */
1115 fprintf(fplog
, " F-max = %12.5e on atom %d\n",
1116 s_min
->fmax
, s_min
->a_fmax
+1);
1117 fprintf(fplog
, " F-Norm = %12.5e\n",
1118 s_min
->fnorm
/sqrtNumAtoms
);
1119 fprintf(fplog
, "\n");
1121 /* Start the loop over CG steps.
1122 * Each successful step is counted, and we continue until
1123 * we either converge or reach the max number of steps.
1126 for (step
= 0; (number_steps
< 0 || (number_steps
>= 0 && step
<= number_steps
)) && !converged
; step
++)
1129 /* start taking steps in a new direction
1130 * First time we enter the routine, beta=0, and the direction is
1131 * simply the negative gradient.
1134 /* Calculate the new direction in p, and the gradient in this direction, gpa */
1139 for (i
= 0; i
< mdatoms
->homenr
; i
++)
1141 if (mdatoms
->cFREEZE
)
1143 gf
= mdatoms
->cFREEZE
[i
];
1145 for (m
= 0; m
< DIM
; m
++)
1147 if (!inputrec
->opts
.nFreeze
[gf
][m
])
1149 p
[i
][m
] = sf
[i
][m
] + beta
*p
[i
][m
];
1150 gpa
-= p
[i
][m
]*sf
[i
][m
];
1151 /* f is negative gradient, thus the sign */
1160 /* Sum the gradient along the line across CPUs */
1163 gmx_sumd(1, &gpa
, cr
);
1166 /* Calculate the norm of the search vector */
1167 get_f_norm_max(cr
, &(inputrec
->opts
), mdatoms
, p
, &pnorm
, NULL
, NULL
);
1169 /* Just in case stepsize reaches zero due to numerical precision... */
1172 stepsize
= inputrec
->em_stepsize
/pnorm
;
1176 * Double check the value of the derivative in the search direction.
1177 * If it is positive it must be due to the old information in the
1178 * CG formula, so just remove that and start over with beta=0.
1179 * This corresponds to a steepest descent step.
1184 step
--; /* Don't count this step since we are restarting */
1185 continue; /* Go back to the beginning of the big for-loop */
1188 /* Calculate minimum allowed stepsize, before the average (norm)
1189 * relative change in coordinate is smaller than precision
1192 for (i
= 0; i
< mdatoms
->homenr
; i
++)
1194 for (m
= 0; m
< DIM
; m
++)
1196 tmp
= fabs(s_min
->s
.x
[i
][m
]);
1205 /* Add up from all CPUs */
1208 gmx_sumd(1, &minstep
, cr
);
1211 minstep
= GMX_REAL_EPS
/sqrt(minstep
/(3*state_global
->natoms
));
1213 if (stepsize
< minstep
)
1219 /* Write coordinates if necessary */
1220 do_x
= do_per_step(step
, inputrec
->nstxout
);
1221 do_f
= do_per_step(step
, inputrec
->nstfout
);
1223 write_em_traj(fplog
, cr
, outf
, do_x
, do_f
, NULL
,
1224 top_global
, inputrec
, step
,
1225 s_min
, state_global
, f_global
);
1227 /* Take a step downhill.
1228 * In theory, we should minimize the function along this direction.
1229 * That is quite possible, but it turns out to take 5-10 function evaluations
1230 * for each line. However, we dont really need to find the exact minimum -
1231 * it is much better to start a new CG step in a modified direction as soon
1232 * as we are close to it. This will save a lot of energy evaluations.
1234 * In practice, we just try to take a single step.
1235 * If it worked (i.e. lowered the energy), we increase the stepsize but
1236 * the continue straight to the next CG step without trying to find any minimum.
1237 * If it didn't work (higher energy), there must be a minimum somewhere between
1238 * the old position and the new one.
1240 * Due to the finite numerical accuracy, it turns out that it is a good idea
1241 * to even accept a SMALL increase in energy, if the derivative is still downhill.
1242 * This leads to lower final energies in the tests I've done. / Erik
1244 s_a
->epot
= s_min
->epot
;
1246 c
= a
+ stepsize
; /* reference position along line is zero */
1248 if (DOMAINDECOMP(cr
) && s_min
->s
.ddp_count
< cr
->dd
->ddp_count
)
1250 em_dd_partition_system(fplog
, step
, cr
, top_global
, inputrec
,
1251 s_min
, top
, mdatoms
, fr
, vsite
, constr
,
1255 /* Take a trial step (new coords in s_c) */
1256 do_em_step(cr
, inputrec
, mdatoms
, fr
->bMolPBC
, s_min
, c
, s_min
->s
.cg_p
, s_c
,
1257 constr
, top
, nrnb
, wcycle
, -1);
1260 /* Calculate energy for the trial step */
1261 evaluate_energy(fplog
, cr
,
1262 top_global
, s_c
, top
,
1263 inputrec
, nrnb
, wcycle
, gstat
,
1264 vsite
, constr
, fcd
, graph
, mdatoms
, fr
,
1265 mu_tot
, enerd
, vir
, pres
, -1, FALSE
);
1267 /* Calc derivative along line */
1271 for (i
= 0; i
< mdatoms
->homenr
; i
++)
1273 for (m
= 0; m
< DIM
; m
++)
1275 gpc
-= p
[i
][m
]*sf
[i
][m
]; /* f is negative gradient, thus the sign */
1278 /* Sum the gradient along the line across CPUs */
1281 gmx_sumd(1, &gpc
, cr
);
1284 /* This is the max amount of increase in energy we tolerate */
1285 tmp
= sqrt(GMX_REAL_EPS
)*fabs(s_a
->epot
);
1287 /* Accept the step if the energy is lower, or if it is not significantly higher
1288 * and the line derivative is still negative.
1290 if (s_c
->epot
< s_a
->epot
|| (gpc
< 0 && s_c
->epot
< (s_a
->epot
+ tmp
)))
1293 /* Great, we found a better energy. Increase step for next iteration
1294 * if we are still going down, decrease it otherwise
1298 stepsize
*= 1.618034; /* The golden section */
1302 stepsize
*= 0.618034; /* 1/golden section */
1307 /* New energy is the same or higher. We will have to do some work
1308 * to find a smaller value in the interval. Take smaller step next time!
1311 stepsize
*= 0.618034;
1317 /* OK, if we didn't find a lower value we will have to locate one now - there must
1318 * be one in the interval [a=0,c].
1319 * The same thing is valid here, though: Don't spend dozens of iterations to find
1320 * the line minimum. We try to interpolate based on the derivative at the endpoints,
1321 * and only continue until we find a lower value. In most cases this means 1-2 iterations.
1323 * I also have a safeguard for potentially really pathological functions so we never
1324 * take more than 20 steps before we give up ...
1326 * If we already found a lower value we just skip this step and continue to the update.
1334 /* Select a new trial point.
1335 * If the derivatives at points a & c have different sign we interpolate to zero,
1336 * otherwise just do a bisection.
1338 if (gpa
< 0 && gpc
> 0)
1340 b
= a
+ gpa
*(a
-c
)/(gpc
-gpa
);
1347 /* safeguard if interpolation close to machine accuracy causes errors:
1348 * never go outside the interval
1350 if (b
<= a
|| b
>= c
)
1355 if (DOMAINDECOMP(cr
) && s_min
->s
.ddp_count
!= cr
->dd
->ddp_count
)
1357 /* Reload the old state */
1358 em_dd_partition_system(fplog
, -1, cr
, top_global
, inputrec
,
1359 s_min
, top
, mdatoms
, fr
, vsite
, constr
,
1363 /* Take a trial step to this new point - new coords in s_b */
1364 do_em_step(cr
, inputrec
, mdatoms
, fr
->bMolPBC
, s_min
, b
, s_min
->s
.cg_p
, s_b
,
1365 constr
, top
, nrnb
, wcycle
, -1);
1368 /* Calculate energy for the trial step */
1369 evaluate_energy(fplog
, cr
,
1370 top_global
, s_b
, top
,
1371 inputrec
, nrnb
, wcycle
, gstat
,
1372 vsite
, constr
, fcd
, graph
, mdatoms
, fr
,
1373 mu_tot
, enerd
, vir
, pres
, -1, FALSE
);
1375 /* p does not change within a step, but since the domain decomposition
1376 * might change, we have to use cg_p of s_b here.
1381 for (i
= 0; i
< mdatoms
->homenr
; i
++)
1383 for (m
= 0; m
< DIM
; m
++)
1385 gpb
-= p
[i
][m
]*sf
[i
][m
]; /* f is negative gradient, thus the sign */
1388 /* Sum the gradient along the line across CPUs */
1391 gmx_sumd(1, &gpb
, cr
);
1396 fprintf(debug
, "CGE: EpotA %f EpotB %f EpotC %f gpb %f\n",
1397 s_a
->epot
, s_b
->epot
, s_c
->epot
, gpb
);
1400 epot_repl
= s_b
->epot
;
1402 /* Keep one of the intervals based on the value of the derivative at the new point */
1405 /* Replace c endpoint with b */
1406 swap_em_state(s_b
, s_c
);
1412 /* Replace a endpoint with b */
1413 swap_em_state(s_b
, s_a
);
1419 * Stop search as soon as we find a value smaller than the endpoints.
1420 * Never run more than 20 steps, no matter what.
1424 while ((epot_repl
> s_a
->epot
|| epot_repl
> s_c
->epot
) &&
1427 if (fabs(epot_repl
- s_min
->epot
) < fabs(s_min
->epot
)*GMX_REAL_EPS
||
1430 /* OK. We couldn't find a significantly lower energy.
1431 * If beta==0 this was steepest descent, and then we give up.
1432 * If not, set beta=0 and restart with steepest descent before quitting.
1442 /* Reset memory before giving up */
1448 /* Select min energy state of A & C, put the best in B.
1450 if (s_c
->epot
< s_a
->epot
)
1454 fprintf(debug
, "CGE: C (%f) is lower than A (%f), moving C to B\n",
1455 s_c
->epot
, s_a
->epot
);
1457 swap_em_state(s_b
, s_c
);
1464 fprintf(debug
, "CGE: A (%f) is lower than C (%f), moving A to B\n",
1465 s_a
->epot
, s_c
->epot
);
1467 swap_em_state(s_b
, s_a
);
1476 fprintf(debug
, "CGE: Found a lower energy %f, moving C to B\n",
1479 swap_em_state(s_b
, s_c
);
1483 /* new search direction */
1484 /* beta = 0 means forget all memory and restart with steepest descents. */
1485 if (nstcg
&& ((step
% nstcg
) == 0))
1491 /* s_min->fnorm cannot be zero, because then we would have converged
1495 /* Polak-Ribiere update.
1496 * Change to fnorm2/fnorm2_old for Fletcher-Reeves
1498 beta
= pr_beta(cr
, &inputrec
->opts
, mdatoms
, top_global
, s_min
, s_b
);
1500 /* Limit beta to prevent oscillations */
1501 if (fabs(beta
) > 5.0)
1507 /* update positions */
1508 swap_em_state(s_min
, s_b
);
1511 /* Print it if necessary */
1516 double sqrtNumAtoms
= sqrt(static_cast<double>(state_global
->natoms
));
1517 fprintf(stderr
, "\rStep %d, Epot=%12.6e, Fnorm=%9.3e, Fmax=%9.3e (atom %d)\n",
1518 step
, s_min
->epot
, s_min
->fnorm
/sqrtNumAtoms
,
1519 s_min
->fmax
, s_min
->a_fmax
+1);
1521 /* Store the new (lower) energies */
1522 upd_mdebin(mdebin
, FALSE
, FALSE
, (double)step
,
1523 mdatoms
->tmass
, enerd
, &s_min
->s
, inputrec
->fepvals
, inputrec
->expandedvals
, s_min
->s
.box
,
1524 NULL
, NULL
, vir
, pres
, NULL
, mu_tot
, constr
);
1526 do_log
= do_per_step(step
, inputrec
->nstlog
);
1527 do_ene
= do_per_step(step
, inputrec
->nstenergy
);
1529 /* Prepare IMD energy record, if bIMD is TRUE. */
1530 IMD_fill_energy_record(inputrec
->bIMD
, inputrec
->imd
, enerd
, step
, TRUE
);
1534 print_ebin_header(fplog
, step
, step
);
1536 print_ebin(mdoutf_get_fp_ene(outf
), do_ene
, FALSE
, FALSE
,
1537 do_log
? fplog
: NULL
, step
, step
, eprNORMAL
,
1538 mdebin
, fcd
, &(top_global
->groups
), &(inputrec
->opts
));
1541 /* Send energies and positions to the IMD client if bIMD is TRUE. */
1542 if (do_IMD(inputrec
->bIMD
, step
, cr
, TRUE
, state_global
->box
, state_global
->x
, inputrec
, 0, wcycle
) && MASTER(cr
))
1544 IMD_send_positions(inputrec
->imd
);
1547 /* Stop when the maximum force lies below tolerance.
1548 * If we have reached machine precision, converged is already set to true.
1550 converged
= converged
|| (s_min
->fmax
< inputrec
->em_tol
);
1552 } /* End of the loop */
1554 /* IMD cleanup, if bIMD is TRUE. */
1555 IMD_finalize(inputrec
->bIMD
, inputrec
->imd
);
1559 step
--; /* we never took that last step in this case */
1562 if (s_min
->fmax
> inputrec
->em_tol
)
1566 warn_step(stderr
, inputrec
->em_tol
, step
-1 == number_steps
, FALSE
);
1567 warn_step(fplog
, inputrec
->em_tol
, step
-1 == number_steps
, FALSE
);
1574 /* If we printed energy and/or logfile last step (which was the last step)
1575 * we don't have to do it again, but otherwise print the final values.
1579 /* Write final value to log since we didn't do anything the last step */
1580 print_ebin_header(fplog
, step
, step
);
1582 if (!do_ene
|| !do_log
)
1584 /* Write final energy file entries */
1585 print_ebin(mdoutf_get_fp_ene(outf
), !do_ene
, FALSE
, FALSE
,
1586 !do_log
? fplog
: NULL
, step
, step
, eprNORMAL
,
1587 mdebin
, fcd
, &(top_global
->groups
), &(inputrec
->opts
));
1591 /* Print some stuff... */
1594 fprintf(stderr
, "\nwriting lowest energy coordinates.\n");
1598 * For accurate normal mode calculation it is imperative that we
1599 * store the last conformation into the full precision binary trajectory.
1601 * However, we should only do it if we did NOT already write this step
1602 * above (which we did if do_x or do_f was true).
1604 do_x
= !do_per_step(step
, inputrec
->nstxout
);
1605 do_f
= (inputrec
->nstfout
> 0 && !do_per_step(step
, inputrec
->nstfout
));
1607 write_em_traj(fplog
, cr
, outf
, do_x
, do_f
, ftp2fn(efSTO
, nfile
, fnm
),
1608 top_global
, inputrec
, step
,
1609 s_min
, state_global
, f_global
);
1614 double sqrtNumAtoms
= sqrt(static_cast<double>(state_global
->natoms
));
1615 fnormn
= s_min
->fnorm
/sqrtNumAtoms
;
1616 print_converged(stderr
, CG
, inputrec
->em_tol
, step
, converged
, number_steps
,
1617 s_min
->epot
, s_min
->fmax
, s_min
->a_fmax
, fnormn
);
1618 print_converged(fplog
, CG
, inputrec
->em_tol
, step
, converged
, number_steps
,
1619 s_min
->epot
, s_min
->fmax
, s_min
->a_fmax
, fnormn
);
1621 fprintf(fplog
, "\nPerformed %d energy evaluations in total.\n", neval
);
1624 finish_em(cr
, outf
, walltime_accounting
, wcycle
);
1626 /* To print the actual number of steps we needed somewhere */
1627 walltime_accounting_set_nsteps_done(walltime_accounting
, step
);
1630 } /* That's all folks */
1633 /*! \brief Do L-BFGS conjugate gradients minimization
1634 \copydoc integrator_t(FILE *fplog, t_commrec *cr,
1635 int nfile, const t_filenm fnm[],
1636 const gmx_output_env_t *oenv, gmx_bool bVerbose,
1638 gmx_vsite_t *vsite, gmx_constr_t constr,
1640 t_inputrec *inputrec,
1641 gmx_mtop_t *top_global, t_fcdata *fcd,
1642 t_state *state_global,
1644 t_nrnb *nrnb, gmx_wallcycle_t wcycle,
1647 int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
1648 gmx_membed_t *membed,
1649 real cpt_period, real max_hours,
1651 unsigned long Flags,
1652 gmx_walltime_accounting_t walltime_accounting)
1654 double do_lbfgs(FILE *fplog
, t_commrec
*cr
,
1655 int nfile
, const t_filenm fnm
[],
1656 const gmx_output_env_t gmx_unused
*oenv
, gmx_bool bVerbose
,
1657 int gmx_unused nstglobalcomm
,
1658 gmx_vsite_t
*vsite
, gmx_constr_t constr
,
1659 int gmx_unused stepout
,
1660 t_inputrec
*inputrec
,
1661 gmx_mtop_t
*top_global
, t_fcdata
*fcd
,
1662 t_state
*state_global
,
1664 t_nrnb
*nrnb
, gmx_wallcycle_t wcycle
,
1665 gmx_edsam_t gmx_unused ed
,
1667 int gmx_unused repl_ex_nst
, int gmx_unused repl_ex_nex
, int gmx_unused repl_ex_seed
,
1668 gmx_membed_t
* /*membed*/,
1669 real gmx_unused cpt_period
, real gmx_unused max_hours
,
1671 unsigned long gmx_unused Flags
,
1672 gmx_walltime_accounting_t walltime_accounting
)
1674 static const char *LBFGS
= "Low-Memory BFGS Minimizer";
1676 gmx_localtop_t
*top
;
1677 gmx_enerdata_t
*enerd
;
1679 gmx_global_stat_t gstat
;
1682 int ncorr
, nmaxcorr
, point
, cp
, neval
, nminstep
;
1683 double stepsize
, step_taken
, gpa
, gpb
, gpc
, tmp
, minstep
;
1684 real
*rho
, *alpha
, *ff
, *xx
, *p
, *s
, *lastx
, *lastf
, **dx
, **dg
;
1685 real
*xa
, *xb
, *xc
, *fa
, *fb
, *fc
, *xtmp
, *ftmp
;
1686 real a
, b
, c
, maxdelta
, delta
;
1687 real diag
, Epot0
, Epot
, EpotA
, EpotB
, EpotC
;
1688 real dgdx
, dgdg
, sq
, yr
, beta
;
1693 gmx_bool do_log
, do_ene
, do_x
, do_f
, foundlower
, *frozen
;
1695 int start
, end
, number_steps
;
1697 int i
, k
, m
, n
, nfmax
, gf
, step
;
1702 gmx_fatal(FARGS
, "Cannot do parallel L-BFGS Minimization - yet.\n");
1707 gmx_fatal(FARGS
, "The combination of constraints and L-BFGS minimization is not implemented. Either do not use constraints, or use another minimizer (e.g. steepest descent).");
1710 n
= 3*state_global
->natoms
;
1711 nmaxcorr
= inputrec
->nbfgscorr
;
1713 /* Allocate memory */
1714 /* Use pointers to real so we dont have to loop over both atoms and
1715 * dimensions all the time...
1716 * x/f are allocated as rvec *, so make new x0/f0 pointers-to-real
1717 * that point to the same memory.
1730 snew(rho
, nmaxcorr
);
1731 snew(alpha
, nmaxcorr
);
1734 for (i
= 0; i
< nmaxcorr
; i
++)
1740 for (i
= 0; i
< nmaxcorr
; i
++)
1749 init_em(fplog
, LBFGS
, cr
, inputrec
,
1750 state_global
, top_global
, &ems
, &top
, &f
, &f_global
,
1751 nrnb
, mu_tot
, fr
, &enerd
, &graph
, mdatoms
, &gstat
, vsite
, constr
,
1752 nfile
, fnm
, &outf
, &mdebin
, imdport
, Flags
, wcycle
);
1753 /* Do_lbfgs is not completely updated like do_steep and do_cg,
1754 * so we free some memory again.
1759 xx
= (real
*)state_global
->x
;
1763 end
= mdatoms
->homenr
;
1765 /* Print to log file */
1766 print_em_start(fplog
, cr
, walltime_accounting
, wcycle
, LBFGS
);
1768 do_log
= do_ene
= do_x
= do_f
= TRUE
;
1770 /* Max number of steps */
1771 number_steps
= inputrec
->nsteps
;
1773 /* Create a 3*natoms index to tell whether each degree of freedom is frozen */
1775 for (i
= start
; i
< end
; i
++)
1777 if (mdatoms
->cFREEZE
)
1779 gf
= mdatoms
->cFREEZE
[i
];
1781 for (m
= 0; m
< DIM
; m
++)
1783 frozen
[3*i
+m
] = inputrec
->opts
.nFreeze
[gf
][m
];
1788 sp_header(stderr
, LBFGS
, inputrec
->em_tol
, number_steps
);
1792 sp_header(fplog
, LBFGS
, inputrec
->em_tol
, number_steps
);
1797 construct_vsites(vsite
, state_global
->x
, 1, NULL
,
1798 top
->idef
.iparams
, top
->idef
.il
,
1799 fr
->ePBC
, fr
->bMolPBC
, cr
, state_global
->box
);
1802 /* Call the force routine and some auxiliary (neighboursearching etc.) */
1803 /* do_force always puts the charge groups in the box and shifts again
1804 * We do not unshift, so molecules are always whole
1807 ems
.s
.x
= state_global
->x
;
1809 evaluate_energy(fplog
, cr
,
1810 top_global
, &ems
, top
,
1811 inputrec
, nrnb
, wcycle
, gstat
,
1812 vsite
, constr
, fcd
, graph
, mdatoms
, fr
,
1813 mu_tot
, enerd
, vir
, pres
, -1, TRUE
);
1818 /* Copy stuff to the energy bin for easy printing etc. */
1819 upd_mdebin(mdebin
, FALSE
, FALSE
, (double)step
,
1820 mdatoms
->tmass
, enerd
, state_global
, inputrec
->fepvals
, inputrec
->expandedvals
, state_global
->box
,
1821 NULL
, NULL
, vir
, pres
, NULL
, mu_tot
, constr
);
1823 print_ebin_header(fplog
, step
, step
);
1824 print_ebin(mdoutf_get_fp_ene(outf
), TRUE
, FALSE
, FALSE
, fplog
, step
, step
, eprNORMAL
,
1825 mdebin
, fcd
, &(top_global
->groups
), &(inputrec
->opts
));
1829 /* This is the starting energy */
1830 Epot
= enerd
->term
[F_EPOT
];
1836 /* Set the initial step.
1837 * since it will be multiplied by the non-normalized search direction
1838 * vector (force vector the first time), we scale it by the
1839 * norm of the force.
1844 double sqrtNumAtoms
= sqrt(static_cast<double>(state_global
->natoms
));
1845 fprintf(stderr
, "Using %d BFGS correction steps.\n\n", nmaxcorr
);
1846 fprintf(stderr
, " F-max = %12.5e on atom %d\n", fmax
, nfmax
+1);
1847 fprintf(stderr
, " F-Norm = %12.5e\n", fnorm
/sqrtNumAtoms
);
1848 fprintf(stderr
, "\n");
1849 /* and copy to the log file too... */
1850 fprintf(fplog
, "Using %d BFGS correction steps.\n\n", nmaxcorr
);
1851 fprintf(fplog
, " F-max = %12.5e on atom %d\n", fmax
, nfmax
+1);
1852 fprintf(fplog
, " F-Norm = %12.5e\n", fnorm
/sqrtNumAtoms
);
1853 fprintf(fplog
, "\n");
1856 // Point is an index to the memory of search directions, where 0 is the first one.
1859 // Set initial search direction to the force (-gradient), or 0 for frozen particles.
1860 for (i
= 0; i
< n
; i
++)
1864 dx
[point
][i
] = ff
[i
]; /* Initial search direction */
1872 // Stepsize will be modified during the search, and actually it is not critical
1873 // (the main efficiency in the algorithm comes from changing directions), but
1874 // we still need an initial value, so estimate it as the inverse of the norm
1875 // so we take small steps where the potential fluctuates a lot.
1876 stepsize
= 1.0/fnorm
;
1878 /* Start the loop over BFGS steps.
1879 * Each successful step is counted, and we continue until
1880 * we either converge or reach the max number of steps.
1885 /* Set the gradient from the force */
1887 for (step
= 0; (number_steps
< 0 || (number_steps
>= 0 && step
<= number_steps
)) && !converged
; step
++)
1890 /* Write coordinates if necessary */
1891 do_x
= do_per_step(step
, inputrec
->nstxout
);
1892 do_f
= do_per_step(step
, inputrec
->nstfout
);
1897 mdof_flags
|= MDOF_X
;
1902 mdof_flags
|= MDOF_F
;
1907 mdof_flags
|= MDOF_IMD
;
1910 mdoutf_write_to_trajectory_files(fplog
, cr
, outf
, mdof_flags
,
1911 top_global
, step
, (real
)step
, state_global
, state_global
, f
, f
);
1913 /* Do the linesearching in the direction dx[point][0..(n-1)] */
1915 /* make s a pointer to current search direction - point=0 first time we get here */
1918 // calculate line gradient in position A
1919 for (gpa
= 0, i
= 0; i
< n
; i
++)
1924 /* Calculate minimum allowed stepsize along the line, before the average (norm)
1925 * relative change in coordinate is smaller than precision
1927 for (minstep
= 0, i
= 0; i
< n
; i
++)
1937 minstep
= GMX_REAL_EPS
/sqrt(minstep
/n
);
1939 if (stepsize
< minstep
)
1945 // Before taking any steps along the line, store the old position
1946 for (i
= 0; i
< n
; i
++)
1953 for (i
= 0; i
< n
; i
++)
1958 /* Take a step downhill.
1959 * In theory, we should find the actual minimum of the function in this
1960 * direction, somewhere along the line.
1961 * That is quite possible, but it turns out to take 5-10 function evaluations
1962 * for each line. However, we dont really need to find the exact minimum -
1963 * it is much better to start a new BFGS step in a modified direction as soon
1964 * as we are close to it. This will save a lot of energy evaluations.
1966 * In practice, we just try to take a single step.
1967 * If it worked (i.e. lowered the energy), we increase the stepsize but
1968 * continue straight to the next BFGS step without trying to find any minimum,
1969 * i.e. we change the search direction too. If the line was smooth, it is
1970 * likely we are in a smooth region, and then it makes sense to take longer
1971 * steps in the modified search direction too.
1973 * If it didn't work (higher energy), there must be a minimum somewhere between
1974 * the old position and the new one. Then we need to start by finding a lower
1975 * value before we change search direction. Since the energy was apparently
1976 * quite rough, we need to decrease the step size.
1978 * Due to the finite numerical accuracy, it turns out that it is a good idea
1979 * to accept a SMALL increase in energy, if the derivative is still downhill.
1980 * This leads to lower final energies in the tests I've done. / Erik
1983 // State "A" is the first position along the line.
1984 // reference position along line is initially zero
1988 // Check stepsize first. We do not allow displacements
1989 // larger than emstep.
1993 // Pick a new position C by adding stepsize to A.
1996 // Calculate what the largest change in any individual coordinate
1997 // would be (translation along line * gradient along line)
1999 for (i
= 0; i
< n
; i
++)
2002 if (delta
> maxdelta
)
2007 // If any displacement is larger than the stepsize limit, reduce the step
2008 if (maxdelta
> inputrec
->em_stepsize
)
2013 while (maxdelta
> inputrec
->em_stepsize
);
2015 // Take a trial step and move the coordinate array xc[] to position C
2016 for (i
= 0; i
< n
; i
++)
2018 xc
[i
] = lastx
[i
] + c
*s
[i
];
2022 // Calculate energy for the trial step in position C
2023 ems
.s
.x
= (rvec
*)xc
;
2025 evaluate_energy(fplog
, cr
,
2026 top_global
, &ems
, top
,
2027 inputrec
, nrnb
, wcycle
, gstat
,
2028 vsite
, constr
, fcd
, graph
, mdatoms
, fr
,
2029 mu_tot
, enerd
, vir
, pres
, step
, FALSE
);
2032 // Calc line gradient in position C
2033 for (gpc
= 0, i
= 0; i
< n
; i
++)
2035 gpc
-= s
[i
]*fc
[i
]; /* f is negative gradient, thus the sign */
2037 /* Sum the gradient along the line across CPUs */
2040 gmx_sumd(1, &gpc
, cr
);
2043 // This is the max amount of increase in energy we tolerate.
2044 // By allowing VERY small changes (close to numerical precision) we
2045 // frequently find even better (lower) final energies.
2046 tmp
= sqrt(GMX_REAL_EPS
)*fabs(EpotA
);
2048 // Accept the step if the energy is lower in the new position C (compared to A),
2049 // or if it is not significantly higher and the line derivative is still negative.
2050 if (EpotC
< EpotA
|| (gpc
< 0 && EpotC
< (EpotA
+tmp
)))
2052 // Great, we found a better energy. We no longer try to alter the
2053 // stepsize, but simply accept this new better position. The we select a new
2054 // search direction instead, which will be much more efficient than continuing
2055 // to take smaller steps along a line. Set fnorm based on the new C position,
2056 // which will be used to update the stepsize to 1/fnorm further down.
2062 // If we got here, the energy is NOT lower in point C, i.e. it will be the same
2063 // or higher than in point A. In this case it is pointless to move to point C,
2064 // so we will have to do more iterations along the same line to find a smaller
2065 // value in the interval [A=0.0,C].
2066 // Here, A is still 0.0, but that will change when we do a search in the interval
2067 // [0.0,C] below. That search we will do by interpolation or bisection rather
2068 // than with the stepsize, so no need to modify it. For the next search direction
2069 // it will be reset to 1/fnorm anyway.
2075 // OK, if we didn't find a lower value we will have to locate one now - there must
2076 // be one in the interval [a,c].
2077 // The same thing is valid here, though: Don't spend dozens of iterations to find
2078 // the line minimum. We try to interpolate based on the derivative at the endpoints,
2079 // and only continue until we find a lower value. In most cases this means 1-2 iterations.
2080 // I also have a safeguard for potentially really pathological functions so we never
2081 // take more than 20 steps before we give up.
2082 // If we already found a lower value we just skip this step and continue to the update.
2086 // Select a new trial point B in the interval [A,C].
2087 // If the derivatives at points a & c have different sign we interpolate to zero,
2088 // otherwise just do a bisection since there might be multiple minima/maxima
2089 // inside the interval.
2090 if (gpa
< 0 && gpc
> 0)
2092 b
= a
+ gpa
*(a
-c
)/(gpc
-gpa
);
2099 /* safeguard if interpolation close to machine accuracy causes errors:
2100 * never go outside the interval
2102 if (b
<= a
|| b
>= c
)
2107 // Take a trial step to point B
2108 for (i
= 0; i
< n
; i
++)
2110 xb
[i
] = lastx
[i
] + b
*s
[i
];
2114 // Calculate energy for the trial step in point B
2115 ems
.s
.x
= (rvec
*)xb
;
2117 evaluate_energy(fplog
, cr
,
2118 top_global
, &ems
, top
,
2119 inputrec
, nrnb
, wcycle
, gstat
,
2120 vsite
, constr
, fcd
, graph
, mdatoms
, fr
,
2121 mu_tot
, enerd
, vir
, pres
, step
, FALSE
);
2125 // Calculate gradient in point B
2126 for (gpb
= 0, i
= 0; i
< n
; i
++)
2128 gpb
-= s
[i
]*fb
[i
]; /* f is negative gradient, thus the sign */
2131 /* Sum the gradient along the line across CPUs */
2134 gmx_sumd(1, &gpb
, cr
);
2137 // Keep one of the intervals [A,B] or [B,C] based on the value of the derivative
2138 // at the new point B, and rename the endpoints of this new interval A and C.
2141 /* Replace c endpoint with b */
2145 /* swap coord pointers b/c */
2155 /* Replace a endpoint with b */
2159 /* swap coord pointers a/b */
2169 * Stop search as soon as we find a value smaller than the endpoints,
2170 * or if the tolerance is below machine precision.
2171 * Never run more than 20 steps, no matter what.
2175 while ((EpotB
> EpotA
|| EpotB
> EpotC
) && (nminstep
< 20));
2177 if (fabs(EpotB
-Epot0
) < GMX_REAL_EPS
|| nminstep
>= 20)
2179 /* OK. We couldn't find a significantly lower energy.
2180 * If ncorr==0 this was steepest descent, and then we give up.
2181 * If not, reset memory to restart as steepest descent before quitting.
2193 /* Search in gradient direction */
2194 for (i
= 0; i
< n
; i
++)
2196 dx
[point
][i
] = ff
[i
];
2198 /* Reset stepsize */
2199 stepsize
= 1.0/fnorm
;
2204 /* Select min energy state of A & C, put the best in xx/ff/Epot
2210 for (i
= 0; i
< n
; i
++)
2221 for (i
= 0; i
< n
; i
++)
2235 for (i
= 0; i
< n
; i
++)
2243 /* Update the memory information, and calculate a new
2244 * approximation of the inverse hessian
2247 /* Have new data in Epot, xx, ff */
2248 if (ncorr
< nmaxcorr
)
2253 for (i
= 0; i
< n
; i
++)
2255 dg
[point
][i
] = lastf
[i
]-ff
[i
];
2256 dx
[point
][i
] *= step_taken
;
2261 for (i
= 0; i
< n
; i
++)
2263 dgdg
+= dg
[point
][i
]*dg
[point
][i
];
2264 dgdx
+= dg
[point
][i
]*dx
[point
][i
];
2269 rho
[point
] = 1.0/dgdx
;
2272 if (point
>= nmaxcorr
)
2278 for (i
= 0; i
< n
; i
++)
2285 /* Recursive update. First go back over the memory points */
2286 for (k
= 0; k
< ncorr
; k
++)
2295 for (i
= 0; i
< n
; i
++)
2297 sq
+= dx
[cp
][i
]*p
[i
];
2300 alpha
[cp
] = rho
[cp
]*sq
;
2302 for (i
= 0; i
< n
; i
++)
2304 p
[i
] -= alpha
[cp
]*dg
[cp
][i
];
2308 for (i
= 0; i
< n
; i
++)
2313 /* And then go forward again */
2314 for (k
= 0; k
< ncorr
; k
++)
2317 for (i
= 0; i
< n
; i
++)
2319 yr
+= p
[i
]*dg
[cp
][i
];
2323 beta
= alpha
[cp
]-beta
;
2325 for (i
= 0; i
< n
; i
++)
2327 p
[i
] += beta
*dx
[cp
][i
];
2337 for (i
= 0; i
< n
; i
++)
2341 dx
[point
][i
] = p
[i
];
2349 /* Test whether the convergence criterion is met */
2350 get_f_norm_max(cr
, &(inputrec
->opts
), mdatoms
, f
, &fnorm
, &fmax
, &nfmax
);
2352 /* Print it if necessary */
2357 double sqrtNumAtoms
= sqrt(static_cast<double>(state_global
->natoms
));
2358 fprintf(stderr
, "\rStep %d, Epot=%12.6e, Fnorm=%9.3e, Fmax=%9.3e (atom %d)\n",
2359 step
, Epot
, fnorm
/sqrtNumAtoms
, fmax
, nfmax
+1);
2361 /* Store the new (lower) energies */
2362 upd_mdebin(mdebin
, FALSE
, FALSE
, (double)step
,
2363 mdatoms
->tmass
, enerd
, state_global
, inputrec
->fepvals
, inputrec
->expandedvals
, state_global
->box
,
2364 NULL
, NULL
, vir
, pres
, NULL
, mu_tot
, constr
);
2365 do_log
= do_per_step(step
, inputrec
->nstlog
);
2366 do_ene
= do_per_step(step
, inputrec
->nstenergy
);
2369 print_ebin_header(fplog
, step
, step
);
2371 print_ebin(mdoutf_get_fp_ene(outf
), do_ene
, FALSE
, FALSE
,
2372 do_log
? fplog
: NULL
, step
, step
, eprNORMAL
,
2373 mdebin
, fcd
, &(top_global
->groups
), &(inputrec
->opts
));
2376 /* Send x and E to IMD client, if bIMD is TRUE. */
2377 if (do_IMD(inputrec
->bIMD
, step
, cr
, TRUE
, state_global
->box
, state_global
->x
, inputrec
, 0, wcycle
) && MASTER(cr
))
2379 IMD_send_positions(inputrec
->imd
);
2382 // Reset stepsize in we are doing more iterations
2383 stepsize
= 1.0/fnorm
;
2385 /* Stop when the maximum force lies below tolerance.
2386 * If we have reached machine precision, converged is already set to true.
2388 converged
= converged
|| (fmax
< inputrec
->em_tol
);
2390 } /* End of the loop */
2392 /* IMD cleanup, if bIMD is TRUE. */
2393 IMD_finalize(inputrec
->bIMD
, inputrec
->imd
);
2397 step
--; /* we never took that last step in this case */
2400 if (fmax
> inputrec
->em_tol
)
2404 warn_step(stderr
, inputrec
->em_tol
, step
-1 == number_steps
, FALSE
);
2405 warn_step(fplog
, inputrec
->em_tol
, step
-1 == number_steps
, FALSE
);
2410 /* If we printed energy and/or logfile last step (which was the last step)
2411 * we don't have to do it again, but otherwise print the final values.
2413 if (!do_log
) /* Write final value to log since we didn't do anythin last step */
2415 print_ebin_header(fplog
, step
, step
);
2417 if (!do_ene
|| !do_log
) /* Write final energy file entries */
2419 print_ebin(mdoutf_get_fp_ene(outf
), !do_ene
, FALSE
, FALSE
,
2420 !do_log
? fplog
: NULL
, step
, step
, eprNORMAL
,
2421 mdebin
, fcd
, &(top_global
->groups
), &(inputrec
->opts
));
2424 /* Print some stuff... */
2427 fprintf(stderr
, "\nwriting lowest energy coordinates.\n");
2431 * For accurate normal mode calculation it is imperative that we
2432 * store the last conformation into the full precision binary trajectory.
2434 * However, we should only do it if we did NOT already write this step
2435 * above (which we did if do_x or do_f was true).
2437 do_x
= !do_per_step(step
, inputrec
->nstxout
);
2438 do_f
= !do_per_step(step
, inputrec
->nstfout
);
2439 write_em_traj(fplog
, cr
, outf
, do_x
, do_f
, ftp2fn(efSTO
, nfile
, fnm
),
2440 top_global
, inputrec
, step
,
2441 &ems
, state_global
, f
);
2445 double sqrtNumAtoms
= sqrt(static_cast<double>(state_global
->natoms
));
2446 print_converged(stderr
, LBFGS
, inputrec
->em_tol
, step
, converged
,
2447 number_steps
, Epot
, fmax
, nfmax
, fnorm
/sqrtNumAtoms
);
2448 print_converged(fplog
, LBFGS
, inputrec
->em_tol
, step
, converged
,
2449 number_steps
, Epot
, fmax
, nfmax
, fnorm
/sqrtNumAtoms
);
2451 fprintf(fplog
, "\nPerformed %d energy evaluations in total.\n", neval
);
2454 finish_em(cr
, outf
, walltime_accounting
, wcycle
);
2456 /* To print the actual number of steps we needed somewhere */
2457 walltime_accounting_set_nsteps_done(walltime_accounting
, step
);
2460 } /* That's all folks */
2462 /*! \brief Do steepest descents minimization
2463 \copydoc integrator_t(FILE *fplog, t_commrec *cr,
2464 int nfile, const t_filenm fnm[],
2465 const gmx_output_env_t *oenv, gmx_bool bVerbose,
2467 gmx_vsite_t *vsite, gmx_constr_t constr,
2469 t_inputrec *inputrec,
2470 gmx_mtop_t *top_global, t_fcdata *fcd,
2471 t_state *state_global,
2473 t_nrnb *nrnb, gmx_wallcycle_t wcycle,
2476 int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
2477 gmx_membed_t *membed,
2478 real cpt_period, real max_hours,
2480 unsigned long Flags,
2481 gmx_walltime_accounting_t walltime_accounting)
2483 double do_steep(FILE *fplog
, t_commrec
*cr
,
2484 int nfile
, const t_filenm fnm
[],
2485 const gmx_output_env_t gmx_unused
*oenv
, gmx_bool bVerbose
,
2486 int gmx_unused nstglobalcomm
,
2487 gmx_vsite_t
*vsite
, gmx_constr_t constr
,
2488 int gmx_unused stepout
,
2489 t_inputrec
*inputrec
,
2490 gmx_mtop_t
*top_global
, t_fcdata
*fcd
,
2491 t_state
*state_global
,
2493 t_nrnb
*nrnb
, gmx_wallcycle_t wcycle
,
2494 gmx_edsam_t gmx_unused ed
,
2496 int gmx_unused repl_ex_nst
, int gmx_unused repl_ex_nex
, int gmx_unused repl_ex_seed
,
2497 gmx_membed_t
* /*membed*/,
2498 real gmx_unused cpt_period
, real gmx_unused max_hours
,
2500 unsigned long gmx_unused Flags
,
2501 gmx_walltime_accounting_t walltime_accounting
)
2503 const char *SD
= "Steepest Descents";
2504 em_state_t
*s_min
, *s_try
;
2506 gmx_localtop_t
*top
;
2507 gmx_enerdata_t
*enerd
;
2509 gmx_global_stat_t gstat
;
2515 gmx_bool bDone
, bAbort
, do_x
, do_f
;
2520 int steps_accepted
= 0;
2522 s_min
= init_em_state();
2523 s_try
= init_em_state();
2525 /* Init em and store the local state in s_try */
2526 init_em(fplog
, SD
, cr
, inputrec
,
2527 state_global
, top_global
, s_try
, &top
, &f
, &f_global
,
2528 nrnb
, mu_tot
, fr
, &enerd
, &graph
, mdatoms
, &gstat
, vsite
, constr
,
2529 nfile
, fnm
, &outf
, &mdebin
, imdport
, Flags
, wcycle
);
2531 /* Print to log file */
2532 print_em_start(fplog
, cr
, walltime_accounting
, wcycle
, SD
);
2534 /* Set variables for stepsize (in nm). This is the largest
2535 * step that we are going to make in any direction.
2537 ustep
= inputrec
->em_stepsize
;
2540 /* Max number of steps */
2541 nsteps
= inputrec
->nsteps
;
2545 /* Print to the screen */
2546 sp_header(stderr
, SD
, inputrec
->em_tol
, nsteps
);
2550 sp_header(fplog
, SD
, inputrec
->em_tol
, nsteps
);
2553 /**** HERE STARTS THE LOOP ****
2554 * count is the counter for the number of steps
2555 * bDone will be TRUE when the minimization has converged
2556 * bAbort will be TRUE when nsteps steps have been performed or when
2557 * the stepsize becomes smaller than is reasonable for machine precision
2562 while (!bDone
&& !bAbort
)
2564 bAbort
= (nsteps
>= 0) && (count
== nsteps
);
2566 /* set new coordinates, except for first step */
2569 do_em_step(cr
, inputrec
, mdatoms
, fr
->bMolPBC
,
2570 s_min
, stepsize
, s_min
->f
, s_try
,
2571 constr
, top
, nrnb
, wcycle
, count
);
2574 evaluate_energy(fplog
, cr
,
2575 top_global
, s_try
, top
,
2576 inputrec
, nrnb
, wcycle
, gstat
,
2577 vsite
, constr
, fcd
, graph
, mdatoms
, fr
,
2578 mu_tot
, enerd
, vir
, pres
, count
, count
== 0);
2582 print_ebin_header(fplog
, count
, count
);
2587 s_min
->epot
= s_try
->epot
;
2590 /* Print it if necessary */
2595 fprintf(stderr
, "Step=%5d, Dmax= %6.1e nm, Epot= %12.5e Fmax= %11.5e, atom= %d%c",
2596 count
, ustep
, s_try
->epot
, s_try
->fmax
, s_try
->a_fmax
+1,
2597 ( (count
== 0) || (s_try
->epot
< s_min
->epot
) ) ? '\n' : '\r');
2600 if ( (count
== 0) || (s_try
->epot
< s_min
->epot
) )
2602 /* Store the new (lower) energies */
2603 upd_mdebin(mdebin
, FALSE
, FALSE
, (double)count
,
2604 mdatoms
->tmass
, enerd
, &s_try
->s
, inputrec
->fepvals
, inputrec
->expandedvals
,
2605 s_try
->s
.box
, NULL
, NULL
, vir
, pres
, NULL
, mu_tot
, constr
);
2607 /* Prepare IMD energy record, if bIMD is TRUE. */
2608 IMD_fill_energy_record(inputrec
->bIMD
, inputrec
->imd
, enerd
, count
, TRUE
);
2610 print_ebin(mdoutf_get_fp_ene(outf
), TRUE
,
2611 do_per_step(steps_accepted
, inputrec
->nstdisreout
),
2612 do_per_step(steps_accepted
, inputrec
->nstorireout
),
2613 fplog
, count
, count
, eprNORMAL
,
2614 mdebin
, fcd
, &(top_global
->groups
), &(inputrec
->opts
));
2619 /* Now if the new energy is smaller than the previous...
2620 * or if this is the first step!
2621 * or if we did random steps!
2624 if ( (count
== 0) || (s_try
->epot
< s_min
->epot
) )
2628 /* Test whether the convergence criterion is met... */
2629 bDone
= (s_try
->fmax
< inputrec
->em_tol
);
2631 /* Copy the arrays for force, positions and energy */
2632 /* The 'Min' array always holds the coords and forces of the minimal
2634 swap_em_state(s_min
, s_try
);
2640 /* Write to trn, if necessary */
2641 do_x
= do_per_step(steps_accepted
, inputrec
->nstxout
);
2642 do_f
= do_per_step(steps_accepted
, inputrec
->nstfout
);
2643 write_em_traj(fplog
, cr
, outf
, do_x
, do_f
, NULL
,
2644 top_global
, inputrec
, count
,
2645 s_min
, state_global
, f_global
);
2649 /* If energy is not smaller make the step smaller... */
2652 if (DOMAINDECOMP(cr
) && s_min
->s
.ddp_count
!= cr
->dd
->ddp_count
)
2654 /* Reload the old state */
2655 em_dd_partition_system(fplog
, count
, cr
, top_global
, inputrec
,
2656 s_min
, top
, mdatoms
, fr
, vsite
, constr
,
2661 /* Determine new step */
2662 stepsize
= ustep
/s_min
->fmax
;
2664 /* Check if stepsize is too small, with 1 nm as a characteristic length */
2666 if (count
== nsteps
|| ustep
< 1e-12)
2668 if (count
== nsteps
|| ustep
< 1e-6)
2673 warn_step(stderr
, inputrec
->em_tol
, count
== nsteps
, constr
!= NULL
);
2674 warn_step(fplog
, inputrec
->em_tol
, count
== nsteps
, constr
!= NULL
);
2679 /* Send IMD energies and positions, if bIMD is TRUE. */
2680 if (do_IMD(inputrec
->bIMD
, count
, cr
, TRUE
, state_global
->box
, state_global
->x
, inputrec
, 0, wcycle
) && MASTER(cr
))
2682 IMD_send_positions(inputrec
->imd
);
2686 } /* End of the loop */
2688 /* IMD cleanup, if bIMD is TRUE. */
2689 IMD_finalize(inputrec
->bIMD
, inputrec
->imd
);
2691 /* Print some data... */
2694 fprintf(stderr
, "\nwriting lowest energy coordinates.\n");
2696 write_em_traj(fplog
, cr
, outf
, TRUE
, inputrec
->nstfout
, ftp2fn(efSTO
, nfile
, fnm
),
2697 top_global
, inputrec
, count
,
2698 s_min
, state_global
, f_global
);
2702 double sqrtNumAtoms
= sqrt(static_cast<double>(state_global
->natoms
));
2703 fnormn
= s_min
->fnorm
/sqrtNumAtoms
;
2705 print_converged(stderr
, SD
, inputrec
->em_tol
, count
, bDone
, nsteps
,
2706 s_min
->epot
, s_min
->fmax
, s_min
->a_fmax
, fnormn
);
2707 print_converged(fplog
, SD
, inputrec
->em_tol
, count
, bDone
, nsteps
,
2708 s_min
->epot
, s_min
->fmax
, s_min
->a_fmax
, fnormn
);
2711 finish_em(cr
, outf
, walltime_accounting
, wcycle
);
2713 /* To print the actual number of steps we needed somewhere */
2714 inputrec
->nsteps
= count
;
2716 walltime_accounting_set_nsteps_done(walltime_accounting
, count
);
2719 } /* That's all folks */
2721 /*! \brief Do normal modes analysis
2722 \copydoc integrator_t(FILE *fplog, t_commrec *cr,
2723 int nfile, const t_filenm fnm[],
2724 const gmx_output_env_t *oenv, gmx_bool bVerbose,
2726 gmx_vsite_t *vsite, gmx_constr_t constr,
2728 t_inputrec *inputrec,
2729 gmx_mtop_t *top_global, t_fcdata *fcd,
2730 t_state *state_global,
2732 t_nrnb *nrnb, gmx_wallcycle_t wcycle,
2735 int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
2736 gmx_membed_t *membed,
2737 real cpt_period, real max_hours,
2739 unsigned long Flags,
2740 gmx_walltime_accounting_t walltime_accounting)
2742 double do_nm(FILE *fplog
, t_commrec
*cr
,
2743 int nfile
, const t_filenm fnm
[],
2744 const gmx_output_env_t gmx_unused
*oenv
, gmx_bool bVerbose
,
2745 int gmx_unused nstglobalcomm
,
2746 gmx_vsite_t
*vsite
, gmx_constr_t constr
,
2747 int gmx_unused stepout
,
2748 t_inputrec
*inputrec
,
2749 gmx_mtop_t
*top_global
, t_fcdata
*fcd
,
2750 t_state
*state_global
,
2752 t_nrnb
*nrnb
, gmx_wallcycle_t wcycle
,
2753 gmx_edsam_t gmx_unused ed
,
2755 int gmx_unused repl_ex_nst
, int gmx_unused repl_ex_nex
, int gmx_unused repl_ex_seed
,
2756 gmx_membed_t
* /*membed*/,
2757 real gmx_unused cpt_period
, real gmx_unused max_hours
,
2759 unsigned long gmx_unused Flags
,
2760 gmx_walltime_accounting_t walltime_accounting
)
2762 const char *NM
= "Normal Mode Analysis";
2764 int natoms
, atom
, d
;
2767 gmx_localtop_t
*top
;
2768 gmx_enerdata_t
*enerd
;
2770 gmx_global_stat_t gstat
;
2775 gmx_bool bSparse
; /* use sparse matrix storage format */
2777 gmx_sparsematrix_t
* sparse_matrix
= NULL
;
2778 real
* full_matrix
= NULL
;
2779 em_state_t
* state_work
;
2781 /* added with respect to mdrun */
2782 int i
, j
, k
, row
, col
;
2783 real der_range
= 10.0*sqrt(GMX_REAL_EPS
);
2785 bool bIsMaster
= MASTER(cr
);
2789 gmx_fatal(FARGS
, "Constraints present with Normal Mode Analysis, this combination is not supported");
2792 state_work
= init_em_state();
2794 /* Init em and store the local state in state_minimum */
2795 init_em(fplog
, NM
, cr
, inputrec
,
2796 state_global
, top_global
, state_work
, &top
,
2798 nrnb
, mu_tot
, fr
, &enerd
, &graph
, mdatoms
, &gstat
, vsite
, constr
,
2799 nfile
, fnm
, &outf
, NULL
, imdport
, Flags
, wcycle
);
2801 natoms
= top_global
->natoms
;
2809 "NOTE: This version of GROMACS has been compiled in single precision,\n"
2810 " which MIGHT not be accurate enough for normal mode analysis.\n"
2811 " GROMACS now uses sparse matrix storage, so the memory requirements\n"
2812 " are fairly modest even if you recompile in double precision.\n\n");
2816 /* Check if we can/should use sparse storage format.
2818 * Sparse format is only useful when the Hessian itself is sparse, which it
2819 * will be when we use a cutoff.
2820 * For small systems (n<1000) it is easier to always use full matrix format, though.
2822 if (EEL_FULL(fr
->eeltype
) || fr
->rlist
== 0.0)
2824 md_print_info(cr
, fplog
, "Non-cutoff electrostatics used, forcing full Hessian format.\n");
2827 else if (top_global
->natoms
< 1000)
2829 md_print_info(cr
, fplog
, "Small system size (N=%d), using full Hessian format.\n", top_global
->natoms
);
2834 md_print_info(cr
, fplog
, "Using compressed symmetric sparse Hessian format.\n");
2840 sz
= DIM
*top_global
->natoms
;
2842 fprintf(stderr
, "Allocating Hessian memory...\n\n");
2846 sparse_matrix
= gmx_sparsematrix_init(sz
);
2847 sparse_matrix
->compressed_symmetric
= TRUE
;
2851 snew(full_matrix
, sz
*sz
);
2859 /* Write start time and temperature */
2860 print_em_start(fplog
, cr
, walltime_accounting
, wcycle
, NM
);
2862 /* fudge nr of steps to nr of atoms */
2863 inputrec
->nsteps
= natoms
*2;
2867 fprintf(stderr
, "starting normal mode calculation '%s'\n%d steps.\n\n",
2868 *(top_global
->name
), (int)inputrec
->nsteps
);
2871 nnodes
= cr
->nnodes
;
2873 /* Make evaluate_energy do a single node force calculation */
2875 evaluate_energy(fplog
, cr
,
2876 top_global
, state_work
, top
,
2877 inputrec
, nrnb
, wcycle
, gstat
,
2878 vsite
, constr
, fcd
, graph
, mdatoms
, fr
,
2879 mu_tot
, enerd
, vir
, pres
, -1, TRUE
);
2880 cr
->nnodes
= nnodes
;
2882 /* if forces are not small, warn user */
2883 get_state_f_norm_max(cr
, &(inputrec
->opts
), mdatoms
, state_work
);
2885 md_print_info(cr
, fplog
, "Maximum force:%12.5e\n", state_work
->fmax
);
2886 if (state_work
->fmax
> 1.0e-3)
2888 md_print_info(cr
, fplog
,
2889 "The force is probably not small enough to "
2890 "ensure that you are at a minimum.\n"
2891 "Be aware that negative eigenvalues may occur\n"
2892 "when the resulting matrix is diagonalized.\n\n");
2895 /***********************************************************
2897 * Loop over all pairs in matrix
2899 * do_force called twice. Once with positive and
2900 * once with negative displacement
2902 ************************************************************/
2904 /* Steps are divided one by one over the nodes */
2905 for (atom
= cr
->nodeid
; atom
< natoms
; atom
+= nnodes
)
2908 for (d
= 0; d
< DIM
; d
++)
2910 x_min
= state_work
->s
.x
[atom
][d
];
2912 state_work
->s
.x
[atom
][d
] = x_min
- der_range
;
2914 /* Make evaluate_energy do a single node force calculation */
2916 evaluate_energy(fplog
, cr
,
2917 top_global
, state_work
, top
,
2918 inputrec
, nrnb
, wcycle
, gstat
,
2919 vsite
, constr
, fcd
, graph
, mdatoms
, fr
,
2920 mu_tot
, enerd
, vir
, pres
, atom
*2, FALSE
);
2922 for (i
= 0; i
< natoms
; i
++)
2924 copy_rvec(state_work
->f
[i
], fneg
[i
]);
2927 state_work
->s
.x
[atom
][d
] = x_min
+ der_range
;
2929 evaluate_energy(fplog
, cr
,
2930 top_global
, state_work
, top
,
2931 inputrec
, nrnb
, wcycle
, gstat
,
2932 vsite
, constr
, fcd
, graph
, mdatoms
, fr
,
2933 mu_tot
, enerd
, vir
, pres
, atom
*2+1, FALSE
);
2934 cr
->nnodes
= nnodes
;
2936 /* x is restored to original */
2937 state_work
->s
.x
[atom
][d
] = x_min
;
2939 for (j
= 0; j
< natoms
; j
++)
2941 for (k
= 0; (k
< DIM
); k
++)
2944 -(state_work
->f
[j
][k
] - fneg
[j
][k
])/(2*der_range
);
2952 #define mpi_type MPI_DOUBLE
2954 #define mpi_type MPI_FLOAT
2956 MPI_Send(dfdx
[0], natoms
*DIM
, mpi_type
, MASTER(cr
), cr
->nodeid
,
2957 cr
->mpi_comm_mygroup
);
2962 for (node
= 0; (node
< nnodes
&& atom
+node
< natoms
); node
++)
2968 MPI_Recv(dfdx
[0], natoms
*DIM
, mpi_type
, node
, node
,
2969 cr
->mpi_comm_mygroup
, &stat
);
2974 row
= (atom
+ node
)*DIM
+ d
;
2976 for (j
= 0; j
< natoms
; j
++)
2978 for (k
= 0; k
< DIM
; k
++)
2984 if (col
>= row
&& dfdx
[j
][k
] != 0.0)
2986 gmx_sparsematrix_increment_value(sparse_matrix
,
2987 row
, col
, dfdx
[j
][k
]);
2992 full_matrix
[row
*sz
+col
] = dfdx
[j
][k
];
2999 if (bVerbose
&& fplog
)
3004 /* write progress */
3005 if (bIsMaster
&& bVerbose
)
3007 fprintf(stderr
, "\rFinished step %d out of %d",
3008 std::min(atom
+nnodes
, natoms
), natoms
);
3015 fprintf(stderr
, "\n\nWriting Hessian...\n");
3016 gmx_mtxio_write(ftp2fn(efMTX
, nfile
, fnm
), sz
, sz
, full_matrix
, sparse_matrix
);
3019 finish_em(cr
, outf
, walltime_accounting
, wcycle
);
3021 walltime_accounting_set_nsteps_done(walltime_accounting
, natoms
*2);