implemented msreact mode (crest-lab#270)

* implemented msreact mode Signed-off-by: gorges <[email protected]> --------- Signed-off-by: gorges <[email protected]> Signed-off-by: Philipp Pracht <[email protected]> Co-authored-by: Philipp Pracht <[email protected]>
pprcht · Mar 6, 2024 · 7e3b360 · 7e3b360
1 parent 89d60ab
commit 7e3b360
Show file tree

Hide file tree

Showing 16 changed files with 1,973 additions and 528 deletions.
diff --git a/.gitignore b/.gitignore
@@ -1,6 +1,8 @@
 *.o
 *.mod
 *.tgz
+*.i90
+*__genmod.f90
 github_bin/
 build_majestix
 build_commands

diff --git a/docs/man/crest.adoc b/docs/man/crest.adoc
@@ -275,6 +275,47 @@ The user is able to include additional constraints to **ALL** xtb**(1)** calcula
   Sum of population for structures considered in msRRHO average.
   [_default_: **0.9** (= 90%)]
 
+=== options for MSREACT automated mass spectra fragment generator
+
+    *-msreact*::
+      start the msreact mode
+
+    *-msnoattrh*::
+      deactivate attractive potential between hydrogen and LMO centers
+
+    *-msnshifts* _int_::  
+      perform n optimizations with randomly shifted atom postions (default 0)
+
+    *-msnshifts* _int_::  
+      perform n optimizations with randomly shifted atom postions and repulsive potential applied to bonds (default 0)
+
+    *-msnbonds* _int_::  
+      maximum number of bonds between atoms pairs for applying repulsive potential (default 3) 
+
+    *-msmolbar*::
+      sort out topological duplicates by molbar codes (requires  sourced "molbar")
+
+    *-msinchi*::
+      sort out topological duplicates by inchi codes (requires  sourced "obabel")
+
+    *-msnfrag* _int_:: 
+        number of fragments that are printed by msreact (random selection)
+
+    *-msiso*::
+      print only non-dissociated structures (isomers)
+
+    *-msnoiso*::
+      print only dissociated structures
+
+    *-mslargeprint*::
+        do not remove temporary files and MSDIR with constrained optimizations
+
+    *-chrg* _int_:: 
+      set the molecules´ charge
+
+    *-ewin* _float_::       
+      set energy window in for sorting out fragments kcal/mol, [default: 200.0 kcal/mol] 
+
 === Other tools for standalone use
 
 *-zsort*::

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
@@ -27,6 +27,7 @@ add_subdirectory("rigidconf")
 add_subdirectory("discretize")
 add_subdirectory("entropy")
 add_subdirectory("legacy_algos")
+add_subdirectory("msreact")
 
 
 set(dir "${CMAKE_CURRENT_SOURCE_DIR}")
@@ -59,8 +60,6 @@ list(APPEND srcs
   "${dir}/marqfit.f90"
   "${dir}/minitools.f90"
   "${dir}/miscdata.f90"
-  "${dir}/msmod.f90"
-  "${dir}/msreact.f90"
   "${dir}/ncigeo.f90"
   "${dir}/ompmklset.F90"
   "${dir}/printouts.f90"

diff --git a/src/classes.f90 b/src/classes.f90
@@ -471,6 +471,19 @@ module crest_data
     type(lwoniom_input),allocatable :: ONIOM_input
     !================================================!
 
+    !>--- msreact mode settings
+    logical :: msnoiso =.false. ! print only dissociated structures in msreact
+    logical :: msiso =.false. ! only print non-dissociated structures in msreact
+    logical :: msmolbar =.false. ! sort out duplicates by molbar
+    logical :: msinchi =.false. ! sort out duplicates by inchi
+    logical :: mslargeprint=.false. ! dont remove temporary files
+    logical :: msattrh=.true. ! add attractive potential for H-atoms
+    integer :: msnbonds = 3 ! distance of bonds up to nonds are stretched
+    integer :: msnshifts = 0 ! number of random shifts applied to whole mol
+    integer :: msnshifts2 = 0 ! number of random shifts applied to whole mol
+    integer :: msnfrag = 0 !number of fragments that are printed in msreact mode
+    character(len=80) :: msinput = '' !name of input file for msreact
+
     !>--- general logical data
     logical :: allrot = .true.       !> use all rotational constants for check instead of mean?
     logical :: altopt = .false.
@@ -567,7 +580,6 @@ module crest_data
     logical :: water = .false.       !> true if water is used as solvent (only QCG)
     logical :: wallsetup = .false.   !> set up a wall potential?
     logical :: wbotopo = .false.     !> set up topo with WBOs
-
   contains
     procedure :: allocate => allocate_metadyn
     procedure :: deallocate => deallocate_metadyn

diff --git a/src/confparse.f90 b/src/confparse.f90
@@ -305,6 +305,17 @@ subroutine parseflags(env,arg,nra)
   env%solv_file = ''
   env%solu_file = ''
 
+!>--- options for msreact
+  env%msiso = .false.  ! msiso and msnoiso are mutually exclusive !!!
+  env%msnoiso = .false.
+  env%msmolbar = .false.
+  env%mslargeprint=.false. ! dont remove temporary files
+  env%msattrh=.true. ! add attractive potential for H-atoms
+  env%msnbonds = 3 ! distance of bonds up to nonds are stretched
+  env%msnshifts = 0 ! number of random shifts applied to whole mol
+  env%msnshifts2 = 0 ! number of random shifts applied to whole mol
+  env%msnfrag = 0 !number of fragments that are printed in msreact mode
+
 !&>
 
 !=========================================================================================!
@@ -536,10 +547,11 @@ subroutine parseflags(env,arg,nra)
         env%autozsort = .false.
         exit
 
-      case ('-msreact')
+    case ('-msreact')
         env%crestver = crest_msreac
         env%preopt = .false.
         env%presp = .true.
+        env%ewin = 200.0d0 !> 200 kcal for msreact
 
       case ('-splitfile')
         ctmp = trim(arg(i+1))
@@ -1007,6 +1019,43 @@ subroutine parseflags(env,arg,nra)
           env%shake = nint(xx(1))
         end select !> QCG
       end if
+
+!========================================================================================!
+!------- Flags for msreact
+!========================================================================================!
+      if (env%crestver == crest_msreac) then
+        select case (argument) !> msreact
+        case('-msnoiso') !> filter out non fragmentated structures in msreact
+          env%msnoiso=.true.
+        case('-msiso') !> filter out fragmentated structures in msreact
+          env%msiso=.true.
+        case('-msnbonds') ! give number of bonds up to which bias potential is added between atoms default 3
+          call readl(arg(i + 1),xx,j)
+          env%msnbonds = xx(1)
+        case('-msnshifts') ! give number of times atoms are randomly shifted before optimization
+          call readl(arg(i + 1),xx,j)
+          env%msnshifts = xx(1)
+        case('-msnshifts2') ! give number of times atoms are randomly shifted before applying the constrained optimization default 0 
+          call readl(arg(i + 1),xx,j)
+          env%msnshifts2 = xx(1)
+        case('-msnfrag') ! give number of structures that should be generated
+          call readl(arg(i + 1),xx,j)
+          env%msnfrag = xx(1)
+        case('-msmolbar') !> filter out structures with same molbar code in msreact
+          env%msmolbar=.true.
+        case('-msinchi') !> filter out structures with same inchi code in msreact
+            env%msinchi=.true.
+        case('-msnoattrh') !> add attractive potential for H-atoms
+          env%msattrh=.false.
+        case('-mslargeprint') !> additional printouts and keep MSDIR
+          env%mslargeprint=.true.
+        case('-msinput') ! give number of times atoms are randomly shifted before applying the constrained optimization default 0 
+            ctmp = trim(arg(i+1))
+            if (ctmp(1:1) .ne. '-') then
+              env%msinput = trim(ctmp)
+            end if
+        end select !> msreact
+      end if
 !========================================================================================!
 !------- other general flags
 !========================================================================================!

diff --git a/src/cregen.f90 b/src/cregen.f90
@@ -28,6 +28,7 @@
 !>              6 - energy sorting only
 !>              9 - no sorting, only check groups
 !>             12 - no topology check, turn ewin to infty
+!>             13 - no topology check, ewin and rmsd checking (msreact settings)
 !=========================================================================================!
 !=========================================================================================!
 
@@ -364,6 +365,11 @@ subroutine cregen_files(env,fname,oname,cname,simpleset,userinput,iounit)
     oname = "crest_mecp_search.xyz.sorted"
     cname = "crest_ensemble.xyz"
   end if
+  if (simpleset == 13) then !> MSREACT files
+    fname = "crest_unique_products.xyz"
+    oname = "crest_unique_products.sorted"
+    cname = "crest_msreact_products.xyz"
+  end if
   if (simpleset == 15) then !> crossing files
     call checkname_xyz('confcross',fname,oname)
     cname = trim(fname)//'.unique'
@@ -420,6 +426,13 @@ subroutine cregen_prout(env,simpleset,pr1,pr2,pr3,pr4)
     pr4 = .true.
   end if
 
+  if (simpleset == 13) then
+    pr1 = .false.
+    pr2 = .false.
+    pr3 = .false.
+    pr4 = .false.
+  end if
+
   return
 end subroutine cregen_prout
 
@@ -515,6 +528,21 @@ subroutine cregen_director(env,simpleset,checkbroken,sortE,sortRMSD,sortRMSD2, &
     anal = .false.
   end if
 
+  if(simpleset == 13)then  !msreact mode 
+    checkbroken = .false.
+    sorte = .true.
+    sortRMSD = .true.
+    sortRMSD2 = .false.
+    repairord = .false.
+    newfile = .true.
+    conffile = .true.
+    topocheck = .false.
+    checkez = .false.
+    bonusfiles= .false.
+    anal = .false.
+
+  endif    
+
   return
 end subroutine cregen_director
 

diff --git a/src/legacy_algos/protonate.f90 b/src/legacy_algos/protonate.f90
@@ -197,7 +197,7 @@ end subroutine protonate
 !--------------------------------------------------------------------------------------------
 ! A quick single point xtb calculation and calculate LMOs
 !--------------------------------------------------------------------------------------------
-subroutine xtblmo(env)
+subroutine xtblmo(env,print)
   use crest_parameters
   use iomod
   use crest_data
@@ -207,6 +207,7 @@ subroutine xtblmo(env)
   character(len=:),allocatable :: jobcall
   integer :: io
   character(len=*),parameter :: pipe = ' > xtb.out 2>/dev/null'
+  logical, optional :: print ! leave the xtb.out file (e.g. for msreact mode)
 
 !---- setting threads
   if (env%autothreads) then
@@ -228,7 +229,7 @@ subroutine xtblmo(env)
 
 !---- cleanup
   call remove(fname)
-  call remove('xtb.out')
+  if (.not. print) call remove('xtb.out')
   call remove('energy')
   call remove('charges')
   call remove('xtbrestart')

diff --git a/src/meson.build b/src/meson.build
@@ -26,6 +26,7 @@ subdir('rigidconf')
 subdir('discretize')
 subdir('entropy')
 subdir('legacy_algos')
+subdir('msreact')
 
 
 srcs += files(
@@ -56,8 +57,6 @@ srcs += files(
   'marqfit.f90',
   'minitools.f90',
   'miscdata.f90',
-  'msmod.f90',
-  'msreact.f90',
   'ncigeo.f90',
   'ompmklset.F90',
   'printouts.f90',

diff --git a/src/miscdata.f90 b/src/miscdata.f90
@@ -69,6 +69,37 @@ module miscdata
    & * aatoau * 4.0_wp / 3.0_wp
 !&>
 
+!&<
+   ! Radius used in QCxMS (in au)
+   real(wp), parameter :: QCxMS_Rad(118) = aatoau *  [ &
+   & 0.32_wp,0.37_wp, & ! H,He
+   & 1.30_wp,0.99_wp,0.84_wp,0.75_wp,0.71_wp,0.64_wp,0.60_wp,0.62_wp, & ! Li-Ne
+   & 1.60_wp,1.40_wp,1.24_wp,1.14_wp,1.09_wp,1.04_wp,1.00_wp,1.01_wp, & ! Na-Ar
+   & 2.00_wp,1.74_wp, & ! K,Ca
+   &                 1.59_wp,1.48_wp,1.44_wp,1.30_wp,1.29_wp, & ! Sc-
+   &                 1.24_wp,1.18_wp,1.17_wp,1.22_wp,1.20_wp, & ! -Zn
+   &                 1.23_wp,1.20_wp,1.20_wp,1.18_wp,1.17_wp,1.16_wp, & ! Ga-Kr
+   & 2.15_wp,1.90_wp, & ! Rb,Sr
+   &                 1.76_wp,1.64_wp,1.56_wp,1.46_wp,1.38_wp, & ! Y-
+   &                 1.36_wp,1.34_wp,1.30_wp,1.36_wp,1.40_wp, & ! -Cd
+   &                 1.42_wp,1.40_wp,1.40_wp,1.37_wp,1.36_wp,1.36_wp, & ! In-Xe
+   & 2.38_wp,2.06_wp, & ! Cs,Ba
+   &         1.94_wp,1.84_wp,1.90_wp,1.88_wp,1.86_wp,1.85_wp,1.83_wp, & ! La-Eu
+   &         1.82_wp,1.81_wp,1.80_wp,1.79_wp,1.77_wp,1.77_wp,1.78_wp, & ! Gd-Yb
+   &                 1.74_wp,1.64_wp,1.58_wp,1.50_wp,1.41_wp, & ! Lu-
+   &                 1.36_wp,1.32_wp,1.30_wp,1.30_wp,1.32_wp, & ! -Hg
+   &                 1.44_wp,1.45_wp,1.50_wp,1.42_wp,1.48_wp,1.46_wp, & ! Tl-Rn
+   & 2.42_wp,2.11_wp, & ! Fr,Ra
+   &         2.01_wp,1.90_wp,1.84_wp,1.83_wp,1.80_wp,1.80_wp,& ! Ac-Pu
+   ! from covalent 2009 covalent radii, such that it is complete up to 118
+   &                                                         1.49_wp, & ! Am
+   &         1.49_wp,1.51_wp,1.51_wp,1.48_wp,1.50_wp,1.56_wp,1.58_wp, & ! Cm-No
+   &                 1.45_wp,1.41_wp,1.34_wp,1.29_wp,1.27_wp, & ! Lr-
+   &                 1.21_wp,1.16_wp,1.15_wp,1.09_wp,1.22_wp, & ! -Cn
+   &                 1.36_wp,1.43_wp,1.46_wp,1.58_wp,1.48_wp,1.57_wp ] ! Nh-Og   
+!&>
+
+
 !&<
   !> D3 pairwise van-der-Waals radii (only homoatomic pairs present here)
   real(wp), parameter :: vdw_D3(1:94) = aatoau * [&