Add MTS test

.github/workflows/gfortran.yml

@@ -244,8 +244,8 @@ jobs:
       CFLAGS: -O0
       GCC_V: ${{ matrix.gcc_v}}
       CODECOV_NAME: ${{format('{0} GCC-{1} OpenMPI-4.0', github.job, matrix.gcc_v)}}
-      OPENMPI_V: "4.0"
-      OPENMPI_PATCH: "0"
+      OPENMPI_V: "4.1"
+      OPENMPI_PATCH: "2"
 
     steps:
     - uses: actions/checkout@v2

sample_inputs/input.in

@@ -9,52 +9,26 @@ irest=0,		! Should we restart from restart.xyz? (ignoring mini.dat)
 
 mdtype='md',		! 'md' = classical MD
 nwalk=1,		! number of random walkers 
-istage=0,		! staging transformation (1), without staging (0)
 nstep=5000,		! number of steps
 dt=20.,			! timestep [au]
 irandom=111313,  	! random seed
-nabin=50,		! what is the gap between ab initio calculations?(RESPA ONLY)
-iqmmm=0,		! 0 - QM/MM OFF, 1 - QM/MM ON
 
 nwrite=10,		! how often some output should be printed (estimators etc.)
 nwritex=5,		! how often should we print coordinates?
 nrest=10,		! how often we print restart files?
 nwritev=0,		! how often we print velocities?
-
-isbc=1,			! spherical boundary conditions
-!rb_sbc=27.5d0,          ! radius of the cluster for sbc
-kb_sbc=0.2d0,		! force constant for spherical potential [au]
-/
-
-&system
-massnames='X','Y'       ! assign atom named X and Y
-masses=16,13            ! with masses 16 and 13.
-! - Monitor the distance between atoms 1 and 2
-ndist=1,		! number of bonds for which we calculate density
-dist1=1,		! first atoms for binning
-dist2=2,		! second atoms for binning (i.e. bond between natfirst[1] and natsecond[1])
-xmin=0.5,		! minimum values for analyzed bonds 
-xmax=2.5,		! maximum values for analyzed bonds
-! - Monitor an angle between atoms 1, 2 and 3
-nang=1,
-ang1=1,		! angle definition: first atoms
-ang2=1,		! angle definition: second atoms
-ang3=1,		! angle definition: third  atoms
-! - Monitor a dihedral angle between atoms 1, 2, 3 and 4
-ndih=1,			! number of analyzed dihedrals
-dih1=1,		! dihedral def: first atoms
-dih2=2,		! dihedral def: second atoms
-dih3=3,		! dihedral def: third atoms
-dih4=4,		! dihedral def: fourth atoms
+nwritef=0,		! how often we print forces?
 /
 
+! Thermostat options
 &nhcopt
 temp=298.15,		! temperature [K] for Maxwell-Boltzmann sampling and thermostat
 inose=0,		! Thermostating: Nose-Hoover 1, microcanonical 0,GLE 2, Langevin 3
 nchain=4,		! number of nose-hoover chains
 tau0=0.0015,		! relaxation time of NHC thermostat
 /
 
+! Surface hopping options
 &sh
 istate_init=2,		! initial electronic state
 nstate=3,		! number of electronic states
@@ -68,73 +42,30 @@ alpha=0.1		! parameter for decoherence correction
 revmom=0,               ! Do not reverse momentum in case of frustrated hops.
 /
 
+&system
+massnames='X','Y'       ! assign atom named X and Y
+masses=16,13            ! with masses 16 and 13.
+! - Monitor the distance between atoms 1 and 2
+ndist=1,		! number of bonds for which we calculate density
+dist1=1,		! first atoms for binning
+dist2=2,		! second atoms for binning (i.e. bond between natfirst[1] and natsecond[1])
+! - Monitor an angle between atoms 1, 2 and 3
+nang=1,
+ang1=1,		! angle definition: first atoms
+ang2=2,		! angle definition: second atoms
+ang3=3,		! angle definition: third  atoms
+! - Monitor a dihedral angle between atoms 1, 2, 3 and 4
+ndih=1,         ! number of analyzed dihedrals
+dih1=1,		! dihedral def: first atoms
+dih2=2,		! dihedral def: second atoms
+dih3=3,		! dihedral def: third atoms
+dih4=4,		! dihedral def: fourth atoms
+/
+
 ! QMMM input read only if iqmmm=1 in section general
 ! This is not well tested!!!
 &qmmm
 natqm=0,		! number of atoms in QM part,set to 0 for MM calculation
 natmm=27,		! number of atoms in MM part
 attypes='O','H','Ar'	! atom types
-q=-0.834,+0.417,0	! atomic charges (TIP3P WATER)
-rmin=3.5364,0.0,3.822,	! L-J atomic radius [A]
-eps=0.0002423919,0.0,0.00037954	! L-J well depth [au]
 /
-
-&nab
-cutoff=100.0,	! non-bonded cutoff (both LJ and coulomb)
-nsnb=10,	! update non-bonded list every nsnb steps	
-
-ipbc=0,		!not tested: turns on periodic boundary conditions
-alpha_pme=0,	! parameter for ewald(TO DO)
-kappa_pme=0,	! parameter for ewald(TO DO)
-!ips=0,		!deprecated:  isotropic periodic sum 1-both LJ and coul, 2-coul 3- LJ
-/
-
-----------------------END OF INPUT---------------------------------------------
-
-
-
-
-------------------OTHER OPTIONS (if needed, move them to their namelists)------- 
-
-GENERAL
-gamm=100,		! first coefficient for minimization
-gammthr=1e-11,		! lowest gamma for minimization
-idebug=0,		! printing out some debugging info
-anal_ext=0,		! do we call user-defined function analyze_ext?
-nproc=1,		! number of processors, only for PIMD
-icv=0,			! do we compute heat capacities? 
-ihess=0,		! do we calculate hessian for heat capacities?
-enmini=0,		! number of steps for which we do not cumulate energy(PIMD only)
-END GENERAL
-
-SYSTEM
-!---------SHAKE OPTIONS  (see manual for complete description)
-nshake=0,		! Number of bond constrains
-shake_tol=0.0001	! Shake tolerance
-ishake1=2,2		! list of bond constraints
-ishake2=7,8		!
-nmol=1,			! number of "molecules"
-nshakemol=0,            ! number of constraints in molecules
-
-!---ANGLE and Dihedral analysis
-nbin=2000,		! number of bins for densities
-nang=2,			! number of analyzed angles
-ang1=1,		! angle definition: first atoms
-ang2=1,		! angle definition: second atoms
-ang3=1,		! angle definition: third  atoms
-ndih=2,			! number of analyzed dihedrals
-dih1=1,		! dihedral def: first atoms
-dih2=2,		! dihedral def: second atoms
-dih3=3,		! dihedral def: third atoms
-dih4=4,		! dihedral def: fourth atoms
-nbin_ang=700,		! number of bins for angles and dihedrals
-shiftdih=360.0d0 	! 0 for dihs between (-180,180); 360 for (0,360)
-
-
-END SYSTEM
-
-NHCOPT
-ams=0.02,		! Nose-Hoove mass (this is rewritten anyway if ipimd=1)
-imasst=0,		! 1 - Massive thermostating , 0- Global thermostatt
-END NHCOPT
-

sample_inputs/input.in.cmd

@@ -23,4 +23,3 @@ tau0=0.001,		! relaxation time of NHC thermostat
 nrespnose=3,		! number of inner respa steps for NH thermostat
 nyosh=7,		! number of steps in suzuki-yoshida scheme,can be 1,3 or 7
 /
-

sample_inputs/input.in.gle

@@ -14,9 +14,7 @@ nrest=5,		! how often we print restart files?
 nwritev=0,		! how often we print velocities?
 /
 
-
 &nhcopt
 temp=298.15,		! temperature [K] for Maxwell-Boltzmann sampling and thermostat
-inose=2,		! Thermostating: Nose-Hoover 1, microcanonical 0,GLE 2
+inose=2,		! GLE thermostat for QT of PI+GLE
 /
-

sample_inputs/input.in.pigle

@@ -17,9 +17,7 @@ nrest=5,		! how often we print restart files?
 nwritev=0,		! how often we print velocities?
 /
 
-
 &nhcopt
 temp=298.15,		! temperature [K] for Maxwell-Boltzmann sampling and thermostat
-inose=2,		! GLE  thermostat 2, needs GLE-A and GLE-C files
+inose=2,		! GLE thermostat 2, needs GLE-A and GLE-C files
 /
-

sample_inputs/input.in.pimd

@@ -1,25 +1,23 @@
-This is a sample input file for ABIN for GLE simulation (quantum thermostat).
+This is a sample input file for PIMD simulation in ABIN,
+using staging transformation and Nose-Hoover Chain thermostat
 
 &general
-pot='g09'		! where do we obtain forces? options are:g09,orca,tera,turbo,molpro,orca,qchem 
+pot='g09'
 mdtype='pimd',		
-nwalk=20,               ! 4 is usually sufficient for PIGLE simulations
-istage=1,               ! istage must be 1 for PIMD
+nwalk=20,               ! number of path integral beads
+istage=1,               ! turn on staging transformation (this must be ON)
 nstep=50000,		! number of steps
 dt=20.,			! timestep [au]
 irandom=1651563,  	! random seed
-irest=0,		! should we restart from restart.xyz? (ignoring mini.dat)
+irest=0,		! should we restart from restart.xyz?
 
 nwrite=1,		! how often some output should be printed (estimators etc.)
 nwritex=5,		! how often should we print coordinates?
 nrest=5,		! how often we print restart files?
-nwritev=0,		! how often we print velocities?
 /
 
-
 &nhcopt
 temp=298.15,		! temperature [K] for Maxwell-Boltzmann sampling and thermostat
-inose=1,		! NHC  thermostat 
-tau0=0.001,             ! 
+inose=1,		! NHC thermostat 
+tau0=0.001,             ! thermostat relaxation time (in picoseconds)
 /
-

sample_inputs/input.in.sh

@@ -16,7 +16,7 @@ nwritev=0,		! how often we print velocities?
 /
 
 &nhcopt
-inose=0,		! Thermostating: Nose-Hoover 1, microcanonical 0,GLE 2, Langevin 3
+inose=0,		! NVE ensemble
 !temp=0.00,		! Usually, you would take initial velocities from WIGNER or set them to zero
 /
 
@@ -29,4 +29,3 @@ EnergyDifThr=0.50,      ! maximum energy difference between two consecutive step
 EnergyDriftThr=0.50,    ! maximum energy drift from initial total energy
 substep=100,		! number of substeps for solving ESCH
 /
-

sample_inputs/input.in.shake

@@ -1,22 +1,17 @@
-This is a sample input file MD with SHAKE constraints
+This is a sample input file for classical MD with SHAKE constraints
 
 &general
-pot='nab'		!not done yet,options should be g09,orca,tera,turbo,molpro,nab,harm,morse,guillot,2dho
-nwalk=1,		!number of random walkers 
-istage=0,		!staging transformation (1), without staging (0)
-nstep=10,		!
-dt=20.,			!number of steps and timestep
-irandom=131313,  	!random seed
+pot='mmwater'
+mdtype='md'
+nstep=10,
+dt=20.,
+irandom=131313,
 
 nwrite=1,		! how often some output should be printed (estimators etc.)
 nwritex=1,		! how often should we print coordinates?
 nwritev=1,
 nrest=1,		! how often we print restart files?
 irest=0,		! should we restart from restart.xyz? (ignoring mini.dat), NOT working yet!!
-
-isbc=1,			! spherical boundary conditions
-!rb_sbc=10.0d0,              ! radius of the cluster for sbc
-kb_sbc=0.01d0,		! force constant for spherical potential,in A
 /
 
 &system
@@ -31,16 +26,10 @@ nmolt=53,               ! number of separate systems(partitions) to couple NHC t
 natmolt=53*3            ! number of atoms in each partition
 nshakemol=2,2,0,2,2,3*0,  ! number of constraints in each partition
 imasst=0,               ! DO NOT use massive thermostat when using SHAKE
-inose=1,		! Thermostating: Nose-Hoover 1, microcanonical 0
+inose=1,		! Only global Nose-Hoover is available for SHAKE
 nchain=4,		! number of nose-hoover chains
 temp=298.15,		! initial temperature for Maxwell-Boltzmann sampling [au]
 tau0=0.001
 nyosh=7,
 nrespnose=3
 /
-
-&nab
-cutoff=100.0,	! non-bonded cutoff (both LJ and coulomb)
-nsnb=100,	!update non-bonded list every nsnb steps	
-/
-

src/force_abin.F90

@@ -232,7 +232,9 @@ subroutine force_abin(x, y, z, fx, fy, fz, eclas, chpot, walkmax)
          close (unit=MAXUNITS + iw + walkmax, status='delete')
       end if
 
-      if (iqmmm == 1) call oniom(x, y, z, fx, fy, fz, eclas, iw)
+      if (iqmmm == 1) then
+         call oniom(x, y, z, fx, fy, fz, eclas, iw)
+      end if
 
    end do
 !$OMP END DO

src/force_bound.F90

@@ -16,7 +16,7 @@ module mod_sbc
 
    subroutine sbc_init(x, y, z)
       use mod_const, only: ANG
-      use mod_general, only: natom !,nwalk
+      use mod_general, only: natom
       use mod_system, only: am, names
       real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)
       real(DP) :: r, rmax, xcm, ycm, zcm
@@ -75,13 +75,14 @@ subroutine sbc_init(x, y, z)
 
    end subroutine
 
-   subroutine force_sbc(x, y, z, fx, fy, fz)
+   subroutine force_sbc(x, y, z, fx, fy, fz, walkmax)
       use mod_const, only: ANG
-      use mod_general, only: natom, nwalk, it, nwrite
+      use mod_general, only: natom, it, nwrite
       use mod_system, only: names
       use mod_files, only: URADIUS
       real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)
       real(DP), intent(inout) :: fx(:, :), fy(:, :), fz(:, :)
+      integer, intent(in) :: walkmax
       real(DP) :: r, frb, rmax, xcm, ycm, zcm
       integer :: iat, iw
 
@@ -93,8 +94,7 @@ subroutine force_sbc(x, y, z, fx, fy, fz)
 
       ! write(*,*)'COM: ',xcm,ycm,zcm
 
-      ! kb_sbc and rb_sbc must be specified in input.in
-      do iw = 1, nwalk
+      do iw = 1, walkmax
          do iat = 1, natom
             r = (x(iat, iw) - xcm)**2 + (y(iat, iw) - ycm)**2 + (z(iat, iw) - zcm)**2
             r = dsqrt(r)