photox/abin

Multipurpose ab initio MD program.

abin's Issues

Add PLUMED to CI build

A note of how it could be done:

https://docs.github.com/en/free-pro-team@latest/actions/learn-github-actions/essential-features-of-github-actions#adding-scripts-to-your-workflow

i.e. a PLUMED build script can be added to .github/scripts/ so that the build complexity is there and the gfortran.yaml can stay simple.

Automatically upload artifacts for failing tests in GHA

This should help tremendeously to debug if a test fails in Github Actions. See

https://github.com/actions/upload-artifact

Include IntelONE API compiler in Github Actions

For reference, see
https://github.com/grimme-lab/crest/blob/master/.github/workflows/build.yml

Depends on #27

Wrong energy from water potentials for PIMD

Energy must be divided by number of beads, preferably in forces.F90 (as opposed to inside the various force functions).
We will need to fix the PIMD test after the fix.

Create basic unit tests

pFUnit testing framework seems to be maintained and could be nice to use.

https://github.com/Goddard-Fortran-Ecosystem/pFUnit

FRUITPy could be an alternative:
https://github.com/acroucher/FRUITPy

Other schemes for FSSH hopping probability

Take a look at the following article for another approach for determining the surface hopping probability.

https://aip.scitation.org/doi/10.1063/5.0024372

Basile thinks it might be interesting to try implementing this and try it out. It shouldn't be too difficult to implement according to the paper authors, and hopefully the recent refactor might help as well.

Add other surface-hopping scheme.

Implement Landau-Zener transition probability as described here:
http://arxiv.org/abs/1412.3047

Doesn't need NA couplings, but second energy derivative.

Also, implement Domcke algorithm for tunnelling.

Investigate restart with NHC thermostat

It looks like restarting a simulation using Nose-Hoover thermostat does not result in exactly same results. Need to investigate deeper.

MPI interface to terachem

Add MPI interface to TeraChem. Might lead to large speed ups for small molecules.
Partly done, it needs to be generalized for PIMD and polished up.
Currently in branch mpitera.

Stabilize tests

We should modify tests so that each test is only 1 or 2 MD steps (+restart). That should minimize the differences between compilers/machines so we could be more strict when comparing against reference values.

Convert all calls to abinerror() to fatal_error()

Add missing tests

Some major functionality is still not tested at all, here's a TODO list:

• Multiple time-step ab initio MD using two potentials
• GLE / PIGLE (tests are currently disabled)
• PIGLET (disabled)
• classical MD with GLE
• analytical potentials (Morse, Double well)
• numerical 1D potential using cubic splines
• plumed with well-tempered metadynamics
• some pieces of Surface Hopping code
• global NHC thermostat

interface to `xtb`

test PIMD and PIGLE

We should do some more test of GLE and PIGLE methods and compare to i-Py.
What about the estimators?
Why does the primitive energy estimator does not work? Two possibilities:

1. Bug in the code.
2. It's not supposed to work, maybe because of the way it is derived.

Validate global NHC thermostat, make it default for classical MD?

Currently, by default we use massive Nosé-Hoover thermostat, where each atom is thermostatted separately (in each cartesian dimension). This setup is crucial to make PIMD work, as it ensures fast thermalization of all degrees of freedom. But for classical MD, it might be suboptimal as it might be too aggressive and hinder sampling of slower motions of the system such as diffusion / barrier crossing etc. The global NHC is implemented as well, but it has been validated much less and I am not sure if many people actually used it for production simulations. We should validate that it works and perhaps make it a default for classical MD.

In particular, we should investigate how global NHC works with respect to sampling translations and rotations. I think it should conserve translational / rotational momentum, in which case if we remove them at the beginning, the system will not sample them. In that case we might need to change the variable f which determines the number of conserved quantities and influences the computation of kinetic temperature.

See also #129.

Division by zero in Landau-Zener subroutine

When running the new LZ_ENE test locally, I noticed this message:

Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG

To find out more, I added the -ffpe-trap=invalid compiler flag to FFLAGS in make.vars.
After recompilation, the LZ_ENE test failed with this message:

Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:
#0  0x2ac4bbea419f in ???
#1  0x52dab6 in __mod_lz_MOD_lz_hop
	at /home/hollas/ABIN-DEV/landau_zener.F90:173
#2  0x4f48fd in abin_dyn
	at /home/hollas/ABIN-DEV/abin.F90:230
#3  0x4f5cb3 in main
	at /home/hollas/ABIN-DEV/abin.F90:20

We probably need some if statement to guard against dividing by zero in landau_zener.F90:173

Use `--implicit-none` compiler option

This option should be a default in the Makefile or in the configure - > make.vars.
This would mean we will automatically catch any code that uses implicit typing in Github Actions in the future.

Convert system calls to `execute_command_line`

system() is a GNU intrinsic, although it does appear to work for ifort. execute_command_line subroutine comes from Fortran 2008 standard and is hopefully supported by recent compilers.

https://gcc.gnu.org/onlinedocs/gfortran/EXECUTE_005fCOMMAND_005fLINE.html#EXECUTE_005fCOMMAND_005fLINE

Better initialization of momenta in GLE-MD

I suspect we do not initialize the atomic momenta correctly when using the GLE thermostat, leading for example to a very slow convergence towards the final temperature in Quantum Thermostat simulations.

See this line in code:

ABIN/gle.F90

Line 305 in 960bb30

! TODO, actually pass this to initialize momenta

Create and impose some code standard

We should make at least some basic rules and try to stick to them.

For white-space, we could use e.g. this:
https://github.com/pseewald/fprettify

See also: http://fortranwiki.org/fortran/show/Tools

Use 'use mpi' instead of 'include mpif.h'

Now that we have automated MPI tests, this is a safe change to do.

More info at:
Goddard-Fortran-Ecosystem/pFUnit#34

Make FFTW conditional.

We should not require user to install FFTW, most of the functionality can do without it (and it is not even used at this point).

Guard against division by zero `try_hop_rescale_nacme`

An interesting edge case (i.e. a bug :-D) in the FSSH hopping subroutine was found when running SH with small deltaE threshold. In that case, a hop between states could happen right between the steps where the energy separation exceeds the threshold, and NACME become zero (because it was not calculated once the energy difference between the states is bigger than deltaE). This leads to a division by zero since we are trying to scale velocities along the zero NAC vector.

ABIN/src/surfacehop.F90

Line 931 in 63d62af

g_temp = (b_temp - dsqrt(c_temp)) / 2.0D0 / a_temp

One solution would be to pass in the interpolated NACM, which should never be zero, since otherwise the hopping probability is also zero. However, that is slightly inconsistent since we're rescaling velocities at the end of the time step, not interpolated velocities.

Alternatively, we could simply guard against this condition and call isotropic rescaling when this happens.

Neither of these fixes is very satisfactory. The fundamental issue is that we're trying to hop every substep, but then the rescaling a recalculation happens at the end of the big step. To make things really consistent I think we should try to hop only once per big time step and accumulate hopping probability along the small time step. I know that both approaches are discussed in the literature, but I don't know what is the current consensus. We should check what SHARC/Newton-X does, and perhaps other SH packages.

Run ABIN through Valgrind

GCC also offers some sanitization options that we might want to enable in the CI.

https://gcc.gnu.org/onlinedocs/gcc/Instrumentation-Options.html

Broken surface hopping code if pot.ne.molpro

It is now possible to have basically arbitrary "pot" parameter. However, part of the SH code is expecting pot='molpro'. This needs to change.

Interface to DL-FIND library.

Add interface to DL-FIND library so that we could efficiently optimize with any programme. It can also do the CI search.

Bad PRNG in create_trajectories.sh

We cannot use the same PRNG for ABIN and for create_trajectories.sh.
That is STUPID obviously!! It might not be disaster, since the PRNG that we use has history, so the trajectories won't probably have exactly same sequences, but it is still dum. We need to implement some other generator for this task and test it with small crush.

Compilation with optimalization points to use of uninitialized variables

We should check whether the Warnings are valid (at least those in transform.F90 seems to be.

Invalid intent for input parameters in `lz_restin`

I got tipped off by a compiler warning. It looks like the position and velocity arrays passed into lz_restin should be intent(in) instead of intent(out).

ABIN/src/landau_zener.F90

Line 515 in 4ca6004

subroutine lz_restin(fileunit, x, y, z, vx, vy, vz)

Note that this is actually a bug, since variables declared as intent(out) are not guaranteed to have the values that you passed in, so you might be assigning random numbers here.

ABIN/src/landau_zener.F90

Line 542 in 4ca6004

x_prev = x; y_prev = y; z_prev = z; vx_prev = vx; vy_prev = vy; vz_prev = vz;

I don't understand why the compiler proceeds without error here, supposedly the purpose of the intent attribute is to prevent such mistakes. 😦

Investigate optimal value of tau0 for NHC thermostat

The current value of the tau0 is 0.001 ps, which is very very low, likely tailored towards the PIMD application (but even there it might be suboptimal). I think for most other MD programs such as Amber or Gromacs the value is more like 1 ps. But it is true they are tailored more towards simulations on nanoseconds (and beyond) timescales, whereas in out ab initio world we're typically reaching tens of picosecond. Shorter value guarantees that the system is equilibrated quickly towards target tempertature, but it might hinder sampling of slow modes / process such as difussion / barrier crossing. We should think about a good test case to verify this and choose a perhaps more reasonable value as the default, at least for classical MD. The other option would be to have no default and leave the choice to the user, but we should still do some testing and provides guidance.

Plumed 2.8.0 fails metadynamics unit test

https://github.com/PHOTOX/ABIN/runs/5409736811?check_suite_focus=true

Metadynamics forces produce non-negligible numerical difference in the unit test. We need to investigate if this is a possible bug in Plumed, or some change in interface or something else. There were apparently some changes in the Fortran interface in 2.8.0, perhaps some of them were not backwards compatible.

ERROR STOP *** Encountered 1 or more failures/errors during testing. ***

Error termination. Backtrace:
Time:         0.032 seconds
  
Failure
 in: 
test_plumed_suite.test_force_metad
  Location: 
[test_plumed.pf:320]
fx(1) changed
AssertEqual failure:
      Expected: <0.99991611203271580>
        Actual: <0.99991594977761633>
    Difference: <-0.16225509946732330E-006> (greater than tolerance of 0.0000000000000000)

Enable canonical sampling with GLE thermostat

The code is already there (inose==4), we just need to enable it and test it. I am not sure whether this ever worked, so we need to test carefully.

Better Makefile

• Streamline and cleanup Makefile
• Move ABIN source to src/ or SRC/ directory
• Rename all file to use .F90 extension

Routine read_nacmrest potentialy unsafe

Routines for reading/wrinting NACM for surface hopping restart are potentially unsafe. They are not printing/reading the full NAC matrix, but they print based on the tocalc() array.

This is potentially unsafe, as the energy can be slightly different after restart and tocalc array can thus change for edge cases.

We should probably always print the full NAC matrix, even if it might get quite big.

Make SHAKE user friendly

Automatically determine bonds to hydrogens and water
Do some rigorous testing: see p. 416 in Understanding mol. simulations
Leave this probably to the next version, it requires coupling to global NHC, which is not well tested

Use `open(newunit=unit..` instead of random integers

This is a Fortran 2008 that should be available in all modern compilers. I already successfully use it in unittests/test_plumed.pf.
This applies to many random places around the codebase, but also in mod_files routines.
As part of that, we should move mod_files from modules.F90 to its own file.

Terachem-MPI for Surface Hopping interface broken

Unfortunately, I broke the surface hopping interface to TC in #62. We pass coordinates in wrong units.

QMMM via NAB interface.

Add QMMM interface via NAB libraries.
Partly done, but segfaults.
Is it work doing.
It might be better to just use semiempirical methods such as PM7 or DFTB as the "MM" part.

Generate test coverage

Apparently, this is possible by Gfortran, but we will need to write a couple processing scripts and collect the coverage data from all TESTS.

http://www.softeng.rl.ac.uk/st/archive/SoftEng/SESP/html/SoftwareTools/gcov.html

Set default value for Nose-Hoover.

We should do some more testing of the best value of tau0 for NHC thermostat.
Test systems: 1 water, small cluster of waters

Add interface for CASPT2/CASSCF dynamics.

Add MOLPRO interface for surface hopping with CASPT2 PES but with CASSCF NACs.
Might be useful, but the user have to be sure that CASSF and CASPT2 states are identical in nature.

This feature might became obsolete when CASPT2 NAC are published in MOLPRO (apparently they are already implemented, see recent Todd's paper on ethylene).

Adaptive time-step for surface hopping

SH routines needs to be revamped. After that one could implement adaptive time-step, similar in spirit to FMS approach. This should allow for larger time-step while also allowing for better energy conservation.

This could be well tested with new TeraChem-FMS MPI interface, which needs to be updated.
In MOLPRO, one needs to add diabatization to r.molpro to allow for monitoring of CI vectors along the trajectories.

PBC and Ewald

EWALD jak je to s 1-4 interakcemi? parametr scee.
Jak zabranit truncation error? Ask Kolafa for help with that.

U PBC simulace se mirne lisi sily mme a mme2...nevime proc. (nevim, zda je tohle aktualni)

pro pripadne PME je potreba modifikovat nbond_box!
funkce mme by nemela mit erfc
a co mme2? jak ted delame hessian?

http://archive.ambermd.org/201007/0619.html
amber11/AmberTools/test/nab/Run.ips

nsnb and cutoff for NAB not working?

In the test suite - NAB_HESS - results aparently does not depend on the cutoff and nsnb parameters, which is weird. Potential big bug.

Migrate from Codecov BASH uploader

Codecov bash uploader is being deprecated
https://docs.codecov.com/docs/about-the-codecov-bash-uploader

New uploader:
https://docs.codecov.com/docs/codecov-uploader#introduction

Usage on GHA
https://github.com/marketplace/actions/codecov

CI via Github Actions

We should run tests automatically on every commit and PR. This should now be relatively easy to setup with Github Actions.
Ideally, we would also test and build with MPI. We should definitely test with different gfortran versions.

Examples:
actions/runner-images#202
https://github.com/fortran-lang/fpm/pull/37/files

Pass MPI variables into tests

Need to figure out how to get proper mpirun into TEST/REMD/test.sh
Similar problems may be for other tests that include MPI.

Normal modes.

MD in normal mode representations.
Essential for PIGLET method and for constraints within PIMD.
Essentially done, but does not work.
We should try use the analytic propagator of normal modes and not the one in force_quantum.f90.
See PIGLE publications for details.

PILE tau0

We need to set default tau0 value or warn the user if its not set. Otherwise the program runs with NaN positions.

Custom types for holding simulation state

Maybe something like this

  type :: energies
     ! Ab initio electronic energy
     real(DP) :: e_pot
     ! Harmonic energy of the PI bead necklace
     real(DP) :: e_pi
  end type

  type :: forces
     ! Ab initio forces
     real(DP), allocatable, dimension(:, :) :: fxc, fyc, fzc
     ! Harmonic PI forces
     real(DP), allocatable, dimension(:, :) :: fxq, fyq, fzq
     ! Difference forces between full and reference potential
     real(DP), allocatable, dimension(:, :) :: fxdiff, fydiff, fzdiff
  end type

  type :: positions
     real(DP), allocatable, dimension(:, :) ::x, y, z
  end type

  type :: momenta
     real(DP), allocatable, dimension(:, :) ::px, py, pz
  end type

  type :: electronic_structure
     ! Dipoles, Transition Dipoles
     ! NAC couplings?
     real(DP), allocatable, dimension(:) :: dip, tdip
  end type

  type mdstate
     type(energies) :: e
     type(forces) :: f
     type(momenta) :: p
     type(positions) :: q
     type(electronic_structure) :: elstruct
  end type

  type(mdstate) :: state, state_prev

For surfaceHopping, we could then easily hold the whole state from a previous step, and be able to just do

state_prev = prev

or vice versa in case we wanted to implement adaptive time step. See also https://stackoverflow.com/questions/19111471/fortran-derived-type-assignment