Apply an external electric field

[Louis Lehmann]

Hello,

I was wondering if it is possible to optimize

the atomic dipoles under the influence of an external

field. I found the function

get_potential_and_electric_field_on_points

in

https://github.com/paesanilab/MBX/blob/master/plugins/python/mbx/mbx.py

but im not shure if it has the desired effect?

Also i need to provide an electric potential.

I would appreciate any help.

Greetings

Louis

[Marc Riera]

Hello Louis,

It is possible, but that feature is stil under testing. I believe it is working well, but I have not written any piece of documentation for that. If you are interested in helping with the test, I would be happy to have a chat with you and explain in detail how to do it. It is not complicated, but one has to have the ideas and steps very clear to do it properly.

What do you think? If you want, shoot me an email on my gmail account and we can see if what you want to do is actually possible with the current version of MBX.

See you!

[Louis Lehmann]

Hi Marc,

maybe I'm missing something but I think I cannot write you a mail

to your gmail account, as its not fully visible to me.

I could write you to your institutional one.

So what I want to do, is to add a static external field in the dipole-SCF field routine, in post processing.

i.e. I basically want to add a static external field in (45) from the following publication:

https://doi.org/10.1063/1.478797,

or whatever technique you are using.

Then I want to get the induced dipole moment in order to

calculate electrostatic observables, e.g. the effective polarizabilities

Delta P = alpha_eff * Delta E_ext

I'm currently only using the MB-Pol water model.

So I guess my request would be only the basic implementation.

Thanks for your help and greetings.

Louis

[Marc Riera]

Hello Lois,

If you want to add a permanent electric field to the system, I believe it can be done as MBX is right now, as long as you don't need the gradients. The issue with the gradients is that I am not sure is correctly implemented so far.

The only requirement for it to work is that you need an electric field AND an electrostatic potential. If you can get those two, and you let me know if you mind using C++, I can guide you on getting this to work.

Let me know and we will get started!

[Marc Riera]

Hello Louis, 

Just to clarify, we do not NEED to use C++. Python will work too. C++ may make things easier to start with, though.

[Louis Lehmann]

Hi Marc,

that is good news.

I do speak both languages, but I would prefer python, as its faster to code.

as I want to apply homogenous external fields, I can get the electric potential by integrating the external field,

for example if the field is along the x-direction:

V = - x E_ext

However, I do not know, how and where I correctly set the external electric field

and potential.

Thanks for you help and greetings,

Louis

[Marc Riera]

Hello!

Since this still is a bit of an "under development", I wil be more than happy to work on this with you. I am gonna follow up with a detailed explanation shortly.

Thanks!

[Marc Riera]

Hello Louis,

It took a bit of time to test and prepare, but here I am going to show you how can you do this. First, you must load the NRG file into the System class:

-------

#include <cmath>
#include <cassert>

#include <iomanip>
#include <iostream>
#include <fstream>
#include <cstring>
#include <stdexcept>
#include <cstdlib>

#include "io_tools/read_nrg.h"
#include "io_tools/write_nrg.h"

#include "bblock/system.h"

//#define PRINT_GRADS
namespace {

static std::vector<bblock::System> systems;

}  // namespace

int main(int argc, char** argv) {
    // Check number of arguments
    if (argc < 2) {
        std::cerr << "Usage: " << argv[0] << " <input.nrg> [mbx.json]" << std::endl;
        return 0;
    }

    // Load the nrg file
    try {
        std::ifstream ifs(argv[1]);

        if (!ifs) {
            throw std::runtime_error("could not open the NRG file");
        }

        tools::ReadNrg(argv[1], systems);
        ifs.close();

    } catch (const std::exception& e) {
        std::cerr << " ** Error ** : " << e.what() << std::endl;
        return 1;
    }

    // Load JSON file
    std::vector<double> box;

    for (size_t i = 0; i < systems.size(); i++) {
        if (argc > 2) {
            systems[i].SetUpFromJson(argv[2]);
        } else {
            systems[i].SetUpFromJson();
        }
    }

------

This code snippet will read a the NRG file and the JSON file form command line arguments of this executable. At this point, you can set the electric fields. You must have a vector of length N_SITES (System.GetNumSites() will give you that number) with the magnitude of the field at the position of each site, and a vector of length 3*N_SITES with the electric field at each position. Let's call them "phi" and "efield". You can set those in the system with:

System.SetExternalElectrostaticPotentialAndFieldInSites(phi,efield)

And then you can just calculate the energy

System.Energy(True)

And if you want to retrieve the atomic dipoles, just call the function:

System.GetDipoles(mu_perm, mu_ind)

Which will fill the mu_perm with the permanent dipoles (q*r) and mu_ind with the induced dipoles. 

If you want to compare dipoles before and after the field is added:

e_nofield = systems[0].Energy(True);

systems[0].GetDipoles(mu_perm_nofield, mu_ind_nofield);

systems[0].SetExternalElectrostaticPotentialAndFieldInSites(phi,efield)

e_field = systems[0].Energy(True);

systems[0].GetDipoles(mu_perm_field, mu_ind_field);

I hope this helps you a bit. I know that the python interface may be easier to code, but you will have more flexibility with the C++. 

Let me know if there are any other issues or questions!

[Louis Lehmann]

Hello,

that is good news.

Thanks and greetings

Louis