*/
class ElectricFieldDimension
{
- public:
- /*! \brief
- * Adds an option section to specify parameters for this field component.
- */
- void initMdpOptions(IOptionsContainerWithSections *options, const char *sectionName)
+public:
+ /*! \brief
+ * Adds an option section to specify parameters for this field component.
+ */
+ void initMdpOptions(IOptionsContainerWithSections* options, const char* sectionName)
+ {
+ auto section = options->addSection(OptionSection(sectionName));
+ section.addOption(RealOption("E0").store(&a_));
+ section.addOption(RealOption("omega").store(&omega_));
+ section.addOption(RealOption("t0").store(&t0_));
+ section.addOption(RealOption("sigma").store(&sigma_));
+ }
+ /*! \brief
+ * Creates mdp parameters for this field component.
+ */
+ void buildMdpOutput(KeyValueTreeObjectBuilder* builder, const std::string& name) const
+ {
+ builder->addUniformArray<real>("electric-field-" + name, { a_, omega_, t0_, sigma_ });
+ }
+
+ /*! \brief Evaluates this field component at given time.
+ *
+ * \param[in] t The time to evualate at
+ * \return The electric field
+ */
+ real evaluate(real t) const
+ {
+ if (sigma_ > 0)
{
- auto section = options->addSection(OptionSection(sectionName));
- section.addOption(RealOption("E0").store(&a_));
- section.addOption(RealOption("omega").store(&omega_));
- section.addOption(RealOption("t0").store(&t0_));
- section.addOption(RealOption("sigma").store(&sigma_));
+ return a_ * (std::cos(omega_ * (t - t0_)) * std::exp(-square(t - t0_) / (2.0 * square(sigma_))));
}
- /*! \brief
- * Creates mdp parameters for this field component.
- */
- void buildMdpOutput(KeyValueTreeObjectBuilder *builder, const std::string &name) const
+ else
{
- builder->addUniformArray<real>("electric-field-" + name, {a_, omega_, t0_, sigma_});
- }
-
- /*! \brief Evaluates this field component at given time.
- *
- * \param[in] t The time to evualate at
- * \return The electric field
- */
- real evaluate(real t) const
- {
- if (sigma_ > 0)
- {
- return a_ * (std::cos(omega_*(t-t0_))
- * std::exp(-square(t-t0_)/(2.0*square(sigma_))));
- }
- else
- {
- return a_ * std::cos(omega_*t);
- }
+ return a_ * std::cos(omega_ * t);
}
+ }
- //! Return the amplitude
- real a() const { return a_; }
-
- private:
- //! Coeffient (V / nm)
- real a_ = 0;
- //! Frequency (1/ps)
- real omega_ = 0;
- //! Central time point (ps) for pulse
- real t0_ = 0;
- //! Width of pulse (ps, if zero there is no pulse)
- real sigma_ = 0;
+ //! Return the amplitude
+ real a() const { return a_; }
+
+private:
+ //! Coeffient (V / nm)
+ real a_ = 0;
+ //! Frequency (1/ps)
+ real omega_ = 0;
+ //! Central time point (ps) for pulse
+ real t0_ = 0;
+ //! Width of pulse (ps, if zero there is no pulse)
+ real sigma_ = 0;
};
/*! \internal
* The electric field can be pulsed and oscillating, simply
* oscillating, or static, in each of X,Y,Z directions.
*/
-class ElectricField final : public IMDModule,
- public IMdpOptionProvider, public IMDOutputProvider,
- public IForceProvider
+class ElectricField final : public IMDModule, public IMdpOptionProvider, public IMDOutputProvider, public IForceProvider
{
- public:
- ElectricField() : fpField_(nullptr) {}
+public:
+ ElectricField() : fpField_(nullptr) {}
- // From IMDModule
- IMdpOptionProvider *mdpOptionProvider() override { return this; }
- IMDOutputProvider *outputProvider() override { return this; }
- void initForceProviders(ForceProviders *forceProviders) override
+ // From IMDModule
+ IMdpOptionProvider* mdpOptionProvider() override { return this; }
+ IMDOutputProvider* outputProvider() override { return this; }
+ void initForceProviders(ForceProviders* forceProviders) override
+ {
+ if (isActive())
{
- if (isActive())
- {
- forceProviders->addForceProvider(this);
- }
+ forceProviders->addForceProvider(this);
}
+ }
- // From IMdpOptionProvider
- void initMdpTransform(IKeyValueTreeTransformRules *transform) override;
- void initMdpOptions(IOptionsContainerWithSections *options) override;
- void buildMdpOutput(KeyValueTreeObjectBuilder *builder) const override;
-
- // From IMDOutputProvider
- void initOutput(FILE *fplog, int nfile, const t_filenm fnm[],
- bool bAppendFiles, const gmx_output_env_t *oenv) override;
- void finishOutput() override;
-
- // From IForceProvider
- //! \copydoc IForceProvider::calculateForces()
- void calculateForces(const ForceProviderInput &forceProviderInput,
- ForceProviderOutput *forceProviderOutput) override;
-
- private:
- //! Return whether or not to apply a field
- bool isActive() const;
-
- /*! \brief Return the field strength
- *
- * \param[in] dim The spatial direction
- * \param[in] t The time (ps)
- * \return The field strength in V/nm units
- */
- real field(int dim, real t) const;
-
- /*! \brief Print the field components to a file
- *
- * \param[in] t The time
- * Will throw and exit with fatal error if file is not open.
- */
- void printComponents(double t) const;
-
- //! The components of the applied electric field in each coordinate dimension
- ElectricFieldDimension efield_[DIM];
- //! File pointer for electric field
- FILE *fpField_;
+ // From IMdpOptionProvider
+ void initMdpTransform(IKeyValueTreeTransformRules* transform) override;
+ void initMdpOptions(IOptionsContainerWithSections* options) override;
+ void buildMdpOutput(KeyValueTreeObjectBuilder* builder) const override;
+
+ // From IMDOutputProvider
+ void initOutput(FILE* fplog, int nfile, const t_filenm fnm[], bool bAppendFiles, const gmx_output_env_t* oenv) override;
+ void finishOutput() override;
+
+ // From IForceProvider
+ //! \copydoc IForceProvider::calculateForces()
+ void calculateForces(const ForceProviderInput& forceProviderInput,
+ ForceProviderOutput* forceProviderOutput) override;
+
+private:
+ //! Return whether or not to apply a field
+ bool isActive() const;
+
+ /*! \brief Return the field strength
+ *
+ * \param[in] dim The spatial direction
+ * \param[in] t The time (ps)
+ * \return The field strength in V/nm units
+ */
+ real field(int dim, real t) const;
+
+ /*! \brief Print the field components to a file
+ *
+ * \param[in] t The time
+ * Will throw and exit with fatal error if file is not open.
+ */
+ void printComponents(double t) const;
+
+ //! The components of the applied electric field in each coordinate dimension
+ ElectricFieldDimension efield_[DIM];
+ //! File pointer for electric field
+ FILE* fpField_;
};
//! Converts dynamic parameters from new mdp format to (E0, omega, t0, sigma).
-void convertParameters(gmx::KeyValueTreeObjectBuilder *builder,
- const std::string &value)
+void convertParameters(gmx::KeyValueTreeObjectBuilder* builder, const std::string& value)
{
const std::vector<std::string> sxt = splitString(value);
if (sxt.empty())
builder->addValue<real>("sigma", fromString<real>(sxt[3]));
}
-void ElectricField::initMdpTransform(IKeyValueTreeTransformRules *rules)
+void ElectricField::initMdpTransform(IKeyValueTreeTransformRules* rules)
{
- rules->addRule().from<std::string>("/electric-field-x").toObject("/electric-field/x")
- .transformWith(&convertParameters);
- rules->addRule().from<std::string>("/electric-field-y").toObject("/electric-field/y")
- .transformWith(&convertParameters);
- rules->addRule().from<std::string>("/electric-field-z").toObject("/electric-field/z")
- .transformWith(&convertParameters);
+ rules->addRule().from<std::string>("/electric-field-x").toObject("/electric-field/x").transformWith(&convertParameters);
+ rules->addRule().from<std::string>("/electric-field-y").toObject("/electric-field/y").transformWith(&convertParameters);
+ rules->addRule().from<std::string>("/electric-field-z").toObject("/electric-field/z").transformWith(&convertParameters);
}
-void ElectricField::initMdpOptions(IOptionsContainerWithSections *options)
+void ElectricField::initMdpOptions(IOptionsContainerWithSections* options)
{
auto section = options->addSection(OptionSection("electric-field"));
efield_[XX].initMdpOptions(§ion, "x");
efield_[ZZ].initMdpOptions(§ion, "z");
}
-void ElectricField::buildMdpOutput(KeyValueTreeObjectBuilder *builder) const
+void ElectricField::buildMdpOutput(KeyValueTreeObjectBuilder* builder) const
{
std::string comment =
- R"(; Electric fields
+ R"(; Electric fields
; Format for electric-field-x, etc. is: four real variables:
; amplitude (V/nm), frequency omega (1/ps), time for the pulse peak (ps),
; and sigma (ps) width of the pulse. Omega = 0 means static field,
efield_[ZZ].buildMdpOutput(builder, "z");
}
-void ElectricField::initOutput(FILE *fplog, int nfile, const t_filenm fnm[],
- bool bAppendFiles, const gmx_output_env_t *oenv)
+void ElectricField::initOutput(FILE* fplog, int nfile, const t_filenm fnm[], bool bAppendFiles, const gmx_output_env_t* oenv)
{
if (isActive())
{
}
else
{
- fpField_ = xvgropen(opt2fn("-field", nfile, fnm),
- "Applied electric field", "Time (ps)",
- "E (V/nm)", oenv);
+ fpField_ = xvgropen(opt2fn("-field", nfile, fnm), "Applied electric field",
+ "Time (ps)", "E (V/nm)", oenv);
}
}
}
bool ElectricField::isActive() const
{
- return (efield_[XX].a() != 0 ||
- efield_[YY].a() != 0 ||
- efield_[ZZ].a() != 0);
+ return (efield_[XX].a() != 0 || efield_[YY].a() != 0 || efield_[ZZ].a() != 0);
}
void ElectricField::printComponents(double t) const
{
- fprintf(fpField_, "%10g %10g %10g %10g\n", t,
- field(XX, t), field(YY, t), field(ZZ, t));
+ fprintf(fpField_, "%10g %10g %10g %10g\n", t, field(XX, t), field(YY, t), field(ZZ, t));
}
-void ElectricField::calculateForces(const ForceProviderInput &forceProviderInput,
- ForceProviderOutput *forceProviderOutput)
+void ElectricField::calculateForces(const ForceProviderInput& forceProviderInput,
+ ForceProviderOutput* forceProviderOutput)
{
if (isActive())
{
- const t_mdatoms &mdatoms = forceProviderInput.mdatoms_;
+ const t_mdatoms& mdatoms = forceProviderInput.mdatoms_;
const double t = forceProviderInput.t_;
- const t_commrec &cr = forceProviderInput.cr_;
+ const t_commrec& cr = forceProviderInput.cr_;
// NOTE: The non-conservative electric field does not have a virial
ArrayRef<RVec> f = forceProviderOutput->forceWithVirial_.force_;
for (int m = 0; (m < DIM); m++)
{
- const real fieldStrength = FIELDFAC*field(m, t);
+ const real fieldStrength = FIELDFAC * field(m, t);
if (fieldStrength != 0)
{
for (index i = 0; i != ssize(f); ++i)
{
// NOTE: Not correct with perturbed charges
- f[i][m] += mdatoms.chargeA[i]*fieldStrength;
+ f[i][m] += mdatoms.chargeA[i] * fieldStrength;
}
}
}
}
}
-} // namespace
+} // namespace
std::unique_ptr<IMDModule> createElectricFieldModule()
{
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff