Code Monkey home page Code Monkey logo

abin's Introduction

DOI CI codecov

What is ABIN?

ABIN is a program for performing ab initio molecular dynamics. It was designed specifically to deal with quantum nuclear effects. It can do path integral simulations and also utilizes quantum thermostat based on General Langevin equation. It can also simulate non-adiabatic events using Surface-hoping algorithm.

The basic philosophy of ABIN program is simple. While the program itself handles the propagation of the system according to the equations of motion, the forces and energies are taken from an external electronic structure program such as ORCA or TeraChem. The call to the chosen external program is handled via a simple shell script interface. Therefore, writing a new interface is rather straightforward and can be done without any changes to ABIN or the ab initio code.

The code is provided under the GNU General Public License. A full text of the license can be found in the file LICENCE.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

A full documentation can be found in the folder DOC.

Installation

To compile the code, you'll need a Fortran and C++ compiler. GNU compilers are tested the most, GFortran and g++ compiler versions >=7.0. Versions >=5.4 likely work as well, but always run the test suite to verify. Intel compiler (ifort) is supported since version >=2018, including the newly open-sourced versions (Intel OneAPI).

The compilation can be as easy as:

$ ./configure && make

Always test the installation by running the test suite:

$ make test

If you modify the source code and want to recompile, you should always clean up before the recompilation:

$ make clean && make

Running the code

After compilation, the executable binary should be available in bin/abin. You can copy the binary somewhere in your PATH, or add bin/ to your PATH. To execute an MD simulation using an input file input.in and initial XYZ coordinates in file geom.xyz, run:

bin/abin -i input.in -x geom.xyz

Run bin/abin --help to see other options. Example input files for various types of simulations can be found in sample_inputs/.

Optional dependencies

Some functionality relies on external libraries. These are optional, and the code automatically recognizes which feature is supported for a given build. Run configure -h to see all the options and how to configure them.

To install the libraries, you can use the install scripts in dev_scripts/. We use these in our Continuous Integration testing suite on Github using the Ubuntu 18.04 image.

The optional libraries are:

  • MPICH: An MPI implementation used for Replica Exchange MD and MPI interface with TeraChem.
    • If you just need REMD you can also use other MPI libraries such as OpenMPI or IntelMPI.
  • FFTW: Fast Fourier Transform library used for normal mode transformation in Path Integral MD.
  • PLUMED: A collection of very useful tools for free energy calculations (MetaDynamics, Umbrella Sampling etc).
  • TCPB-CPP: [EXPERIMENTAL] TCPB interface to TeraChem

Structure of the repository

Path Description
src/ ABIN source code
sample_inputs Sample input files.
interfaces/ BASH interfaces to common ab initio codes.
utils/ Handy scripts that might be useful in conjuction with the MD code.
unit_tests/ Unit tests; run by make unittest (needs pFUnit library installed)
tests/ End-to-End tests; run by make e2etest
dev_scripts/ Setup for ABIN devs and install scripts for optional libraries.

abin's People

Contributors

danielhollas avatar janosjiri avatar juraskov avatar mbarnfield63 avatar stepan-srsen avatar suchanj avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar

abin's Issues

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

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).

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

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.

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

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.

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.

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.

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.

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.

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.

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).

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.

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.

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

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.

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.

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)

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.

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.

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.

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.

Better Makefile

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

PILE tau0

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

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

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.

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).

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.

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. 😦

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

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.