2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2015,2016, 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.
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.
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.
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.
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.
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.
37 * Declares data structure and utilities for electric fields
39 * \author David van der Spoel <david.vanderspoel@icm.uu.se>
40 * \ingroup module_applied_forces
44 #include "electricfield.h"
48 #include "gromacs/commandline/filenm.h"
49 #include "gromacs/fileio/gmxfio.h"
50 #include "gromacs/fileio/gmxfio-xdr.h"
51 #include "gromacs/fileio/readinp.h"
52 #include "gromacs/fileio/warninp.h"
53 #include "gromacs/fileio/xvgr.h"
54 #include "gromacs/gmxlib/network.h"
55 #include "gromacs/math/units.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/mdtypes/commrec.h"
58 #include "gromacs/mdtypes/forcerec.h"
59 #include "gromacs/mdtypes/inputrec.h"
60 #include "gromacs/mdtypes/mdatom.h"
61 #include "gromacs/options/basicoptions.h"
62 #include "gromacs/options/ioptionscontainerwithsections.h"
63 #include "gromacs/options/optionsection.h"
64 #include "gromacs/utility/compare.h"
65 #include "gromacs/utility/exceptions.h"
66 #include "gromacs/utility/fatalerror.h"
67 #include "gromacs/utility/keyvaluetreebuilder.h"
68 #include "gromacs/utility/keyvaluetreetransform.h"
69 #include "gromacs/utility/pleasecite.h"
70 #include "gromacs/utility/strconvert.h"
71 #include "gromacs/utility/stringutil.h"
72 #include "gromacs/utility/txtdump.h"
81 * \brief Declaration of storage unit for fields
83 class ElectricFieldData
86 ElectricFieldData() : a_(0), omega_(0), t0_(0), sigma_(0)
91 * Adds an option section to specify parameters for this field component.
93 void initMdpOptions(IOptionsContainerWithSections
*options
, const char *sectionName
)
95 auto section
= options
->addSection(OptionSection(sectionName
));
96 section
.addOption(RealOption("E0").store(&a_
));
97 section
.addOption(RealOption("omega").store(&omega_
));
98 section
.addOption(RealOption("t0").store(&t0_
));
99 section
.addOption(RealOption("sigma").store(&sigma_
));
102 /*! \brief Evaluates this field component at given time.
104 * \param[in] t The time to evualate at
105 * \return The electric field
107 real
evaluate(real t
) const
111 return a_
* (std::cos(omega_
*(t
-t0_
))
112 * std::exp(-square(t
-t0_
)/(2.0*square(sigma_
))));
116 return a_
* std::cos(omega_
*t
);
120 /*! \brief Initiate the field values
122 * \param[in] a Amplitude
123 * \param[in] omega Frequency
124 * \param[in] t0 Peak of the pulse
125 * \param[in] sigma Width of the pulse
127 void setField(real a
, real omega
, real t0
, real sigma
)
135 //! Return the amplitude
136 real
a() const { return a_
; }
137 //! Return the frequency
138 real
omega() const { return omega_
; }
139 //! Return the time for the peak of the pulse
140 real
t0() const { return t0_
; }
141 //! Return the width of the pulse (0 means inifinite)
142 real
sigma() const { return sigma_
; }
145 //! Coeffient (V / nm)
149 //! Central time point (ps) for pulse
151 //! Width of pulse (ps, if zero there is no pulse)
156 * \brief Describe time dependent electric field
158 * Class that implements a force to be evaluated in mdrun.
159 * The electric field can be pulsed and oscillating, simply
160 * oscillating, or static, in each of X,Y,Z directions.
162 class ElectricField
: public IInputRecExtension
, public IForceProvider
165 ElectricField() : fpField_(nullptr) {}
167 // From IInputRecExtension
168 virtual void doTpxIO(t_fileio
*fio
, bool bRead
);
169 virtual void initMdpTransform(IKeyValueTreeTransformRules
*transform
);
170 virtual void initMdpOptions(IOptionsContainerWithSections
*options
);
171 virtual void broadCast(const t_commrec
*cr
);
172 virtual void compare(FILE *fp
,
173 const IInputRecExtension
*field2
,
176 virtual void printParameters(FILE *fp
, int indent
);
177 virtual void initOutput(FILE *fplog
, int nfile
, const t_filenm fnm
[],
178 bool bAppendFiles
, const gmx_output_env_t
*oenv
);
179 virtual void finishOutput();
180 virtual void initForcerec(t_forcerec
*fr
);
182 //! \copydoc gmx::IForceProvider::calculateForces
183 virtual void calculateForces(const t_commrec
*cr
,
184 const t_mdatoms
*atoms
,
185 PaddedRVecVector
*force
,
189 //! Return whether or not to apply a field
190 bool isActive() const;
192 /*! \brief Add a component to the electric field
194 * The electric field has three spatial dimensions that are
195 * added to the data structure one at a time.
196 * \param[in] dim Dimension, XX, YY, ZZ (0, 1, 2)
197 * \param[in] a Amplitude of the field in V/nm
198 * \param[in] omega Frequency (1/ps)
199 * \param[in] t0 Time of pulse peak (ps)
200 * \param[in] sigma Width of peak (ps)
202 void setFieldTerm(int dim
, real a
, real omega
, real t0
, real sigma
);
204 /*! \brief Return the field strength
206 * \param[in] dim The spatial direction
207 * \param[in] t The time (ps)
208 * \return The field strength in V/nm units
210 real
field(int dim
, real t
) const;
212 /*! \brief Return amplitude of field
214 * \param[in] dim Direction of the field (XX, YY, ZZ)
215 * \return Amplitude of the field
217 real
a(int dim
) const { return efield_
[dim
].a(); }
218 /*! \brief Return frequency of field (1/ps)
220 * \param[in] dim Direction of the field (XX, YY, ZZ)
221 * \return Frequency of the field
223 real
omega(int dim
) const { return efield_
[dim
].omega(); }
224 /*! \brief Return time of pulse peak
226 * \param[in] dim Direction of the field (XX, YY, ZZ)
227 * \return Time of pulse peak
229 real
t0(int dim
) const { return efield_
[dim
].t0(); }
230 /*! \brief Return width of the pulse
232 * \param[in] dim Direction of the field (XX, YY, ZZ)
233 * \return Width of the pulse
235 real
sigma(int dim
) const { return efield_
[dim
].sigma(); }
237 /*! \brief Print the field components to a file
239 * \param[in] t The time
240 * Will throw and exit with fatal error if file is not open.
242 void printComponents(double t
) const;
244 //! The field strength in each dimension
245 ElectricFieldData efield_
[DIM
];
246 //! File pointer for electric field
250 void ElectricField::doTpxIO(t_fileio
*fio
, bool bRead
)
252 // The content of the tpr file for this feature has
253 // been the same since gromacs 4.0 that was used for
255 for (int j
= 0; (j
< DIM
); j
++)
261 if (omega(j
) != 0 || sigma(j
) != 0 || t0(j
) != 0)
266 gmx_fio_do_int(fio
, n
);
267 gmx_fio_do_int(fio
, nt
);
268 std::vector
<real
> aa
, phi
, at
, phit
;
272 phi
.push_back(t0(j
));
273 at
.push_back(omega(j
));
274 phit
.push_back(sigma(j
));
283 gmx_fio_ndo_real(fio
, aa
.data(), n
);
284 gmx_fio_ndo_real(fio
, phi
.data(), n
);
285 gmx_fio_ndo_real(fio
, at
.data(), nt
);
286 gmx_fio_ndo_real(fio
, phit
.data(), nt
);
289 setFieldTerm(j
, aa
[0], at
[0], phi
[0], phit
[0]);
292 gmx_fatal(FARGS
, "Can not handle tpr files with more than one electric field term per direction.");
298 //! Converts static parameters from mdp format to E0.
299 real
convertStaticParameters(const std::string
&value
)
301 // TODO: Better context for the exceptions here (possibly
302 // also convert them to warning_errors or such).
303 const std::vector
<std::string
> sx
= splitString(value
);
308 const int n
= fromString
<int>(sx
[0]);
315 GMX_THROW(InvalidInputError("Only one electric field term supported for each dimension"));
319 GMX_THROW(InvalidInputError("Expected exactly one electric field amplitude value"));
321 return fromString
<real
>(sx
[1]);
324 //! Converts dynamic parameters from mdp format to (omega, t0, sigma).
325 void convertDynamicParameters(gmx::KeyValueTreeObjectBuilder
*builder
,
326 const std::string
&value
)
328 const std::vector
<std::string
> sxt
= splitString(value
);
333 const int n
= fromString
<int>(sxt
[0]);
339 GMX_THROW(InvalidInputError("Please specify 1 omega 0 for non-pulsed fields"));
341 builder
->addValue
<real
>("omega", fromString
<real
>(sxt
[1]));
346 GMX_THROW(InvalidInputError("Please specify 1 omega 0 t0 0 sigma 0 for pulsed fields"));
348 builder
->addValue
<real
>("omega", fromString
<real
>(sxt
[1]));
349 builder
->addValue
<real
>("t0", fromString
<real
>(sxt
[3]));
350 builder
->addValue
<real
>("sigma", fromString
<real
>(sxt
[5]));
353 GMX_THROW(InvalidInputError("Incomprehensible input for electric field"));
357 void ElectricField::initMdpTransform(IKeyValueTreeTransformRules
*rules
)
359 rules
->addRule().from
<std::string
>("/E-x").to
<real
>("/electric-field/x/E0")
360 .transformWith(&convertStaticParameters
);
361 rules
->addRule().from
<std::string
>("/E-xt").toObject("/electric-field/x")
362 .transformWith(&convertDynamicParameters
);
363 rules
->addRule().from
<std::string
>("/E-y").to
<real
>("/electric-field/y/E0")
364 .transformWith(&convertStaticParameters
);
365 rules
->addRule().from
<std::string
>("/E-yt").toObject("/electric-field/y")
366 .transformWith(&convertDynamicParameters
);
367 rules
->addRule().from
<std::string
>("/E-z").to
<real
>("/electric-field/z/E0")
368 .transformWith(&convertStaticParameters
);
369 rules
->addRule().from
<std::string
>("/E-zt").toObject("/electric-field/z")
370 .transformWith(&convertDynamicParameters
);
373 void ElectricField::initMdpOptions(IOptionsContainerWithSections
*options
)
375 //CTYPE ("Format is E0 (V/nm), omega (1/ps), t0 (ps), sigma (ps) ");
376 auto section
= options
->addSection(OptionSection("electric-field"));
377 efield_
[XX
].initMdpOptions(§ion
, "x");
378 efield_
[YY
].initMdpOptions(§ion
, "y");
379 efield_
[ZZ
].initMdpOptions(§ion
, "z");
382 void ElectricField::broadCast(const t_commrec
*cr
)
384 rvec a1
, omega1
, sigma1
, t01
;
388 // Load the parameters read from tpr into temp vectors
389 for (int m
= 0; m
< DIM
; m
++)
392 omega1
[m
] = omega(m
);
393 sigma1
[m
] = sigma(m
);
397 // Broadcasting the parameters
398 gmx_bcast(DIM
*sizeof(a1
[0]), a1
, cr
);
399 gmx_bcast(DIM
*sizeof(omega1
[0]), omega1
, cr
);
400 gmx_bcast(DIM
*sizeof(t01
[0]), t01
, cr
);
401 gmx_bcast(DIM
*sizeof(sigma1
[0]), sigma1
, cr
);
403 // And storing them locally
406 for (int m
= 0; m
< DIM
; m
++)
408 setFieldTerm(m
, a1
[m
], omega1
[m
], t01
[m
], sigma1
[m
]);
413 void ElectricField::initOutput(FILE *fplog
, int nfile
, const t_filenm fnm
[],
414 bool bAppendFiles
, const gmx_output_env_t
*oenv
)
418 please_cite(fplog
, "Caleman2008a");
420 // Optional outpuf file showing the field, see manual.
421 if (opt2bSet("-field", nfile
, fnm
))
425 fpField_
= gmx_fio_fopen(opt2fn("-field", nfile
, fnm
), "a+");
429 fpField_
= xvgropen(opt2fn("-field", nfile
, fnm
),
430 "Applied electric field", "Time (ps)",
437 void ElectricField::finishOutput()
439 if (fpField_
!= nullptr)
441 // This is opened sometimes with xvgropen, sometimes with
442 // gmx_fio_fopen, so we use the least common denominator for closing.
443 gmx_fio_fclose(fpField_
);
448 void ElectricField::initForcerec(t_forcerec
*fr
)
452 fr
->bF_NoVirSum
= TRUE
;
457 void ElectricField::printParameters(FILE *fp
, int indent
)
459 const char *const dimension
[DIM
] = { "X", "Y", "Z" };
460 indent
= pr_title(fp
, indent
, "ElectricField");
461 for (int m
= 0; m
< DIM
; m
++)
463 pr_indent(fp
, indent
);
464 fprintf(fp
, "-%s E0 = %g omega = %g t0 = %g sigma = %g\n",
465 dimension
[m
], a(m
), omega(m
), t0(m
), sigma(m
));
469 void ElectricField::setFieldTerm(int dim
, real a
, real omega
, real t0
, real sigma
)
471 range_check(dim
, 0, DIM
);
472 efield_
[dim
].setField(a
, omega
, t0
, sigma
);
475 real
ElectricField::field(int dim
, real t
) const
477 return efield_
[dim
].evaluate(t
);
480 bool ElectricField::isActive() const
482 return (efield_
[XX
].a() != 0 ||
483 efield_
[YY
].a() != 0 ||
484 efield_
[ZZ
].a() != 0);
487 void ElectricField::compare(FILE *fp
,
488 const IInputRecExtension
*other
,
492 GMX_ASSERT(dynamic_cast<const ElectricField
*>(other
) != nullptr,
493 "Invalid other type");
494 const ElectricField
*f2
= static_cast<const ElectricField
*>(other
);
495 for (int m
= 0; (m
< DIM
); m
++)
499 sprintf(buf
, "inputrec->field[%d]", m
);
500 cmp_real(fp
, buf
, -1, a(m
), f2
->a(m
), reltol
, abstol
);
501 cmp_real(fp
, buf
, -1, omega(m
), f2
->omega(m
), reltol
, abstol
);
502 cmp_real(fp
, buf
, -1, t0(m
), f2
->t0(m
), reltol
, abstol
);
503 cmp_real(fp
, buf
, -1, sigma(m
), f2
->sigma(m
), reltol
, abstol
);
507 void ElectricField::printComponents(double t
) const
509 fprintf(fpField_
, "%10g %10g %10g %10g\n", t
,
510 field(XX
, t
), field(YY
, t
), field(ZZ
, t
));
513 void ElectricField::calculateForces(const t_commrec
*cr
,
514 const t_mdatoms
*mdatoms
,
515 PaddedRVecVector
*force
,
520 rvec
*f
= as_rvec_array(force
->data());
522 for (int m
= 0; (m
< DIM
); m
++)
524 real Ext
= FIELDFAC
*field(m
, t
);
528 // TODO: Check parallellism
529 for (int i
= 0; i
< mdatoms
->homenr
; ++i
)
531 // NOTE: Not correct with perturbed charges
532 f
[i
][m
] += mdatoms
->chargeA
[i
]*Ext
;
536 if (MASTER(cr
) && fpField_
!= nullptr)
545 std::unique_ptr
<IInputRecExtension
> createElectricFieldModule()
547 return std::unique_ptr
<IInputRecExtension
>(new ElectricField());