7a82459f2de3ef96819db47e6954a429e3888533
[alexxy/gromacs.git] / src / gromacs / modularsimulator / shellfcelement.cpp
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2019, by the GROMACS development team, led by
5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6  * and including many others, as listed in the AUTHORS file in the
7  * top-level source directory and at http://www.gromacs.org.
8  *
9  * GROMACS is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public License
11  * as published by the Free Software Foundation; either version 2.1
12  * of the License, or (at your option) any later version.
13  *
14  * GROMACS is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17  * Lesser General Public License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with GROMACS; if not, see
21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
23  *
24  * If you want to redistribute modifications to GROMACS, please
25  * consider that scientific software is very special. Version
26  * control is crucial - bugs must be traceable. We will be happy to
27  * consider code for inclusion in the official distribution, but
28  * derived work must not be called official GROMACS. Details are found
29  * in the README & COPYING files - if they are missing, get the
30  * official version at http://www.gromacs.org.
31  *
32  * To help us fund GROMACS development, we humbly ask that you cite
33  * the research papers on the package. Check out http://www.gromacs.org.
34  */
35 /*! \libinternal
36  * \brief Defines the shell / flex constraints element for the modular simulator
37  *
38  * \author Pascal Merz <pascal.merz@me.com>
39  * \ingroup module_modularsimulator
40  */
41
42 #include "gmxpre.h"
43
44 #include "shellfcelement.h"
45
46 #include "gromacs/domdec/mdsetup.h"
47 #include "gromacs/math/vectypes.h"
48 #include "gromacs/mdlib/constr.h"
49 #include "gromacs/mdlib/force.h"
50 #include "gromacs/mdlib/force_flags.h"
51 #include "gromacs/mdlib/mdatoms.h"
52 #include "gromacs/mdrun/shellfc.h"
53 #include "gromacs/mdtypes/inputrec.h"
54 #include "gromacs/topology/atoms.h"
55
56 #include "energyelement.h"
57 #include "statepropagatordata.h"
58
59 struct gmx_edsam;
60 struct gmx_enfrot;
61 struct gmx_multisim_t;
62 class history_t;
63 struct t_graph;
64
65 namespace gmx
66 {
67 bool ShellFCElement::doShellsOrFlexConstraints(
68         const gmx_mtop_t *mtop, int nflexcon)
69 {
70     if (nflexcon != 0)
71     {
72         return true;
73     }
74     std::array<int, eptNR> n = countPtypes(nullptr, mtop);
75     return n[eptShell] != 0;
76 }
77
78 ShellFCElement::ShellFCElement(
79         StatePropagatorData *statePropagatorData,
80         EnergyElement       *energyElement,
81         bool                 isVerbose,
82         bool                 isDynamicBox,
83         FILE                *fplog,
84         const t_commrec     *cr,
85         const t_inputrec    *inputrec,
86         const MDAtoms       *mdAtoms,
87         t_nrnb              *nrnb,
88         t_forcerec          *fr,
89         t_fcdata            *fcd,
90         gmx_wallcycle       *wcycle,
91         MdScheduleWorkload  *mdScheduleWork,
92         gmx_vsite_t         *vsite,
93         ImdSession          *imdSession,
94         pull_t              *pull_work,
95         Constraints         *constr,
96         const gmx_mtop_t    *globalTopology) :
97     nextNSStep_(-1),
98     nextEnergyCalculationStep_(-1),
99     nextVirialCalculationStep_(-1),
100     nextFreeEnergyCalculationStep_(-1),
101     statePropagatorData_(statePropagatorData),
102     energyElement_(energyElement),
103     localTopology_(nullptr),
104     isDynamicBox_(isDynamicBox),
105     isVerbose_(isVerbose),
106     nSteps_(0),
107     ddBalanceRegionHandler_(cr),
108     fplog_(fplog),
109     cr_(cr),
110     inputrec_(inputrec),
111     mdAtoms_(mdAtoms),
112     nrnb_(nrnb),
113     wcycle_(wcycle),
114     fr_(fr),
115     vsite_(vsite),
116     imdSession_(imdSession),
117     pull_work_(pull_work),
118     fcd_(fcd),
119     mdScheduleWork_(mdScheduleWork),
120     constr_(constr)
121 {
122     shellfc_ = init_shell_flexcon(
123                 fplog, globalTopology,
124                 constr_ ? constr_->numFlexibleConstraints() : 0,
125                 inputrec->nstcalcenergy, DOMAINDECOMP(cr));
126
127     GMX_ASSERT(shellfc_, "ShellFCElement built, but init_shell_flexcon returned a nullptr");
128
129     if (!DOMAINDECOMP(cr))
130     {
131         // This was done in mdAlgorithmsSetupAtomData(), but shellfc
132         // won't be available outside this element.
133         make_local_shells(cr, mdAtoms->mdatoms(), shellfc_);
134     }
135 }
136
137 void ShellFCElement::scheduleTask(
138         Step step, Time time,
139         const RegisterRunFunctionPtr &registerRunFunction)
140 {
141     unsigned int flags = (
142             GMX_FORCE_STATECHANGED |
143             GMX_FORCE_ALLFORCES |
144             (isDynamicBox_ ? GMX_FORCE_DYNAMICBOX : 0) |
145             (nextVirialCalculationStep_ == step ? GMX_FORCE_VIRIAL : 0) |
146             (nextEnergyCalculationStep_ == step ? GMX_FORCE_ENERGY : 0) |
147             (nextFreeEnergyCalculationStep_ == step ? GMX_FORCE_DHDL : 0));
148
149     (*registerRunFunction)(
150             std::make_unique<SimulatorRunFunction>(
151                     [this, step, time, flags]()
152                     {run(step, time, flags); }));
153     nSteps_++;
154 }
155
156 void ShellFCElement::elementSetup()
157 {
158     GMX_ASSERT(localTopology_, "Setup called before local topology was set.");
159
160 }
161
162 void ShellFCElement::run(Step step, Time time, unsigned int flags)
163 {
164     // Disabled functionality
165     gmx_multisim_t *ms               = nullptr;
166     gmx_enfrot     *enforcedRotation = nullptr;
167     t_graph        *graph            = nullptr;
168
169     auto            x      = statePropagatorData_->positionsView();
170     auto            v      = statePropagatorData_->velocitiesView();
171     auto            forces = statePropagatorData_->forcesView();
172     auto            box    = statePropagatorData_->constBox();
173     auto            lambda = ArrayRef<real>();     // disabled
174     history_t      *hist   = nullptr;              // disabled
175
176     tensor          force_vir = {{0}};
177     relax_shell_flexcon(fplog_, cr_, ms, isVerbose_,
178                         enforcedRotation, step,
179                         inputrec_, imdSession_, pull_work_, step == nextNSStep_,
180                         static_cast<int>(flags), localTopology_,
181                         constr_, energyElement_->enerdata(), fcd_,
182                         statePropagatorData_->localNumAtoms(),
183                         x, v, box, lambda, hist, forces, force_vir,
184                         mdAtoms_->mdatoms(), nrnb_, wcycle_, graph,
185                         shellfc_, fr_, mdScheduleWork_, time,
186                         energyElement_->muTot(), vsite_,
187                         ddBalanceRegionHandler_);
188     energyElement_->addToForceVirial(force_vir, step);
189 }
190
191 void ShellFCElement::elementTeardown()
192 {
193     done_shellfc(fplog_, shellfc_, nSteps_);
194 }
195
196 void ShellFCElement::setTopology(const gmx_localtop_t *top)
197 {
198     localTopology_ = top;
199 }
200
201 SignallerCallbackPtr ShellFCElement::registerNSCallback()
202 {
203     return std::make_unique<SignallerCallback>(
204             [this](Step step, Time gmx_unused time)
205             {this->nextNSStep_ = step; });
206 }
207
208 SignallerCallbackPtr ShellFCElement::
209     registerEnergyCallback(EnergySignallerEvent event)
210 {
211     if (event == EnergySignallerEvent::energyCalculationStep)
212     {
213         return std::make_unique<SignallerCallback>(
214                 [this](Step step, Time)
215                 {nextEnergyCalculationStep_ = step; });
216     }
217     if (event == EnergySignallerEvent::virialCalculationStep)
218     {
219         return std::make_unique<SignallerCallback>(
220                 [this](Step step, Time){nextVirialCalculationStep_ = step; });
221     }
222     if (event == EnergySignallerEvent::freeEnergyCalculationStep)
223     {
224         return std::make_unique<SignallerCallback>(
225                 [this](Step step, Time){nextFreeEnergyCalculationStep_ = step; });
226     }
227     return nullptr;
228 }
229 }  // namespace std