Energies and forces for a homemade code

[Javier Hernández Rojas]

Hello developers!

Is it possible to calculate energies and forces for water clusters as could be done in older versions (MBPOL)? That is, without using any interface such as LAMMPS or i-PI. I have seen that in Fortran it uses the command ‘call get_energy (coord,n_at,Vpot,grads)’. Is it necessary to have a .ngr file prepared together with the ‘json_file’ file before making this call?

Thanks in advance

Javier

[Philip Zhu]

All functions defined in MBX repo: `src/main/bblock/system.cpp` are available to call from homemade cpp codes. Fortrain availability of function calls might be limited, but the basic usage as the example you showed works.

Under MBX repo: `src/tests/single_point.cpp`, you can find cpp examples of making gradient calls.

It is easier to have a .nrg file together with json file, but if you use cpp, it is also possible to load the system manually. You need to create a System object, and create Monomers with appropriate coordinates and names, then call AddMonomer. Some example lines of manually loading systems can be found at MBX repo: `src/main/mb_decomp.cpp` (starting from line 93 https://github.com/paesanilab/MBX/blob/d2fee399283ee5e477098544a212c473a2069807/src/main/mb_decomp.cpp#L92-L93)

[Javier Hernández Rojas]

Thanks for your prompt reply.

My codes are in Fortran, and I have to evaluate the energy and force of the system many times, so the .ngr file would change every step. 

As I understand it, the command ‘call get_energy_g’ would be executed, taking the information from the .ngr and json files 

located in the directory where the code is running. Does the command "call initialize_system ()" generate the .ngr file?

Regards,

Javier

[Philip Zhu]

For the example from MBX repo: `examples/PEFs/001_mbpol/fortran/`, with the Fortran code locates in `examples/PEFs/src/fortran`, I believe you do not need .nrg file. It reads an "input", an "input.xyz" file, and the json file.

initialize_system should just prepare internal data structures of what MBX needs. The "input" file provides monomer topologies information, together with "json" file those two are only needed once during initialize_system, if you are evaluating the same system many times.

"input.xyz" provides coordinates, which you don't have to read from file. You can provide coordinates generated within the fortran code.

get_energy as explained in the comment calculates energy with/without gradient, with/without periodic boundary condition, depending on which function you call. All function calls available to Fortran & python can be found here: `src/bblock/external_call.cpp`.

If you are evaluating pbc, box info is passed as parameter to function call of get_energy_pbc_ and you don't need to worry about updating json file's box information, because json is only used during initialization.

[Javier Hernández Rojas]

Many thanks Philip. It seems that now it works. 

On the other hand, I have tried  to do a minimization using the ‘optimize’ executable 

for a water dimer, and the system is taking a long time to reach the minimum (> 10m ??). 

What could I be doing wrong? I attach the input files.

Best,

Javier

[Philip Zhu]

Your input seems completely fine. With your input, it took me less than 1s to optimize. One possibility the I can think of is: sometimes, the executable tries to use the maximum possible number of OpenMP threads, which might significantly deteriorate the performance. You can check if this is the case by doing `top` while the job is running; and the cpu usage should not exceed 100%. If you see cpu usage > 100%, you can fix this issue by running `export OMP_NUM_THREADS=1` before running the `optimize` executable.

The reason for this slow down, specifically for some older machines, is that the number of virtual cores (e.g. 24) might be two times the number of physical cores (e.g. 12) due to hyper-threading, which don’t always work well. Your system default might try to run with 24 threads on 12 cores, which over stresses the system.

If this is not the reason, it may be more complicated system specific issues.

The following log shows the output when I tried to run your input.

```

$ time ~/software/mbx-dev+lammps/mbpol2023_gcc/MBX/bin/optimize input.nrg mbx.json 
________________________________________________________________
 MIN:          It=    0  nfunc=     1  E=        -0.75961 ( 1.75e+00)
 MIN:          It=    1  nfunc=     8  E=        -1.12166 ( 1.53e+00)
 MIN:          It=    2  nfunc=    13  E=        -1.35944 ( 2.03e+00)
 MIN:          It=    3  nfunc=    18  E=        -1.53257 ( 2.56e+00)
 MIN:          It=    4  nfunc=    23  E=        -1.69986 ( 2.95e+00)
 MIN:          It=    5  nfunc=    27  E=        -1.90540 ( 6.10e+00)
 MIN:          It=    6  nfunc=    31  E=        -2.17602 ( 1.22e+00)
 MIN:          It=    7  nfunc=    36  E=        -2.28716 ( 1.57e+00)
 MIN:          It=    8  nfunc=    40  E=        -2.38655 ( 2.28e+00)
 MIN:          It=    9  nfunc=    44  E=        -2.44157 ( 5.42e-01)
 MIN:          It=   10  nfunc=    48  E=        -2.49142 ( 1.00e+00)
 MIN:          It=   11  nfunc=    54  E=        -2.76282 ( 5.82e+00)
 MIN:          It=   12  nfunc=    59  E=        -3.18755 ( 4.38e+00)
 MIN:          It=   13  nfunc=    63  E=        -3.48236 ( 2.66e+00)
 MIN:          It=   14  nfunc=    67  E=        -3.63823 ( 1.73e+00)
 MIN:          It=   15  nfunc=    74  E=        -3.98382 ( 4.26e+00)
 MIN:          It=   16  nfunc=    78  E=        -4.13441 ( 9.52e-01)
 MIN:          It=   17  nfunc=    89  E=        -4.27101 ( 1.89e+00)
 MIN:          It=   18  nfunc=    95  E=        -4.36980 ( 2.05e+00)
 MIN:          It=   19  nfunc=    99  E=        -4.49724 ( 3.06e+00)
 MIN:          It=   20  nfunc=   103  E=        -4.60751 ( 9.01e-01)
 MIN:          It=   21  nfunc=   110  E=        -4.68956 ( 2.15e+00)
 MIN:          It=   22  nfunc=   114  E=        -4.76499 ( 1.63e+00)
 MIN:          It=   23  nfunc=   118  E=        -4.80494 ( 6.90e-01)
 MIN:          It=   24  nfunc=   122  E=        -4.82983 ( 1.61e+00)
 MIN:          It=   25  nfunc=   129  E=        -4.87948 ( 1.11e+00)
 MIN:          It=   26  nfunc=   134  E=        -4.89654 ( 6.33e-01)
 MIN:          It=   27  nfunc=   138  E=        -4.90602 ( 5.57e-01)
 MIN:          It=   28  nfunc=   142  E=        -4.91140 ( 2.84e-01)
 MIN:          It=   29  nfunc=   151  E=        -4.93454 ( 9.23e-01)
 MIN:          It=   30  nfunc=   155  E=        -4.94614 ( 3.46e-01)
 MIN:          It=   31  nfunc=   162  E=        -4.95409 ( 5.62e-01)
 MIN:          It=   32  nfunc=   166  E=        -4.95603 ( 1.08e-01)
 MIN:          It=   33  nfunc=   172  E=        -4.95905 ( 4.13e-01)
 MIN:          It=   34  nfunc=   176  E=        -4.95981 ( 6.78e-02)
 MIN:          It=   35  nfunc=   185  E=        -4.96199 ( 2.95e-01)
 MIN:          It=   36  nfunc=   191  E=        -4.96297 ( 2.95e-02)
 MIN:          It=   37  nfunc=   199  E=        -4.96321 ( 8.82e-03)
 MIN:          It=   38  nfunc=   207  E=        -4.96323 ( 7.48e-03)
 MIN:          It=   39  nfunc=   213  E=        -4.96324 ( 6.97e-03)
 MIN:          It=   40  nfunc=   218  E=        -4.96324 ( 1.83e-03)
 MIN:          It=   41  nfunc=   228  E=        -4.96324 ( 5.01e-04)
 MIN:          It=   42  nfunc=   234  E=        -4.96324 ( 8.35e-04)
 MIN:          It=   43  nfunc=   242  E=        -4.96324 ( 3.33e-04)
 MIN:          It=   44  nfunc=   249  E=        -4.96324 ( 3.33e-04)
 MIN:          It=   45  nfunc=   254  E=        -4.96324 ( 2.82e-04)
 MIN:          It=   46  nfunc=   263  E=        -4.96324 ( 4.05e-05)
 MIN:          It=   47  nfunc=   269  E=        -4.96324 ( 4.50e-05)
 MIN:          It=   48  nfunc=   276  E=        -4.96324 ( 2.04e-05)
 MIN:          It=   49  nfunc=   282  E=        -4.96324 ( 1.78e-06)
 MIN:          It=   50  nfunc=   287  E=        -4.96324 ( 2.72e-06)
 MIN:          It=   51  nfunc=   296  E=        -4.96324 ( 1.55e-07)
----------------------------------------------------------------
 FIN:             :-)                    E=      -4.96 (  0.000)

Optimization Converged!
E = -0.759606 -> -4.96324, |g| = 3.10509e-07
-4.96324181e+00   ++++!!!

real    0m0.351s
user    0m0.329s
sys     0m0.017s

$ ls -l
total 26
-rw-rw-r--. 1 xuz144 paesanilab 329 Jan 27 14:52 input.nrg
-rw-rw-r--. 1 xuz144 paesanilab 801 Jan 27 14:52 mbx.json
-rw-rw-r--. 1 xuz144 paesanilab 414 Jan 27 14:53 optimized.xyz
-rw-rw-r--. 1 xuz144 paesanilab   0 Jan 27 14:53 unconverged.xyz

[Javier Hernández Rojas]

Hi developers,

I have performed a single_point calculation with two water molecules, one Li+ atom and one Cl- atom. 

The code determines an energy of -78.56 kcal/mol. Is this energy correct? That is, is it due to n-body interactions 

(2B,3B, pol) except for the ion-ion interaction? Thanks in advance. I attach the input files.

Best,

Javier