fixes (start v3.0.1) (crest-lab#286)

CMakeLists.txt

@@ -27,10 +27,9 @@ endif()
 project(
   crest
   LANGUAGES "C" "Fortran"
-  VERSION  3.0 
+  VERSION  3.0.1 
   DESCRIPTION "A tool for the exploration of low-energy chemical space"
 )
-#set(SOVERSION "pre")
 
 # Follow GNU conventions for installing directories
 include(GNUInstallDirs)

config/CMakeLists.txt

@@ -108,6 +108,12 @@ configure_file(
   "${PROJECT_BINARY_DIR}/crest_metadata.fh"
   @ONLY
 )
+# just to be safe, create it also in include/
+configure_file(
+  "${PROJECT_SOURCE_DIR}/assets/template/metadata.f90"
+  "${PROJECT_BINARY_DIR}/include/crest_metadata.fh"
+  @ONLY
+)
 
 #########################################################################################
 #########################################################################################

meson.build

@@ -17,7 +17,7 @@
 project(
   'crest',
   'fortran', 'c',
-  version: '3.0',
+  version: '3.0.1',
   license: 'LGPL-3.0-or-later',
   meson_version: '>=0.63',
   default_options: [

src/algos/hessian_tools.f90

@@ -44,8 +44,11 @@ module hessian_tools
 !=========================================================================================!
 !=========================================================================================!
 
-  !Returns the Frequencies from a Hessian in cm-1
+
   subroutine frequencies(nat,at,xyz,nat3,calc,prj_mw_hess,freq,io)
+!*************************************************
+!* Returns the Frequencies from a Hessian in cm-1
+!*************************************************
     implicit none
 
     integer,intent(in) :: nat
@@ -93,11 +96,12 @@ subroutine frequencies(nat,at,xyz,nat3,calc,prj_mw_hess,freq,io)
 
   end subroutine frequencies
 
+
   subroutine mass_weight_hess(nat,at,nat3,hess)
     use atmasses
     implicit none
 
-    !Mass weighting the Hessian
+    !> Mass weighting the Hessian
     integer,intent(in) :: nat                   !Number of atoms
     integer,intent(in) :: at(nat)               !atomic number of all atoms
 
@@ -128,8 +132,9 @@ subroutine mass_weight_hess(nat,at,nat3,hess)
   end subroutine mass_weight_hess
 
   subroutine dsqtoh(n,a,b)
-
-    !converts upper triangle of a matrix into a vector
+!****************************************************
+!* converts upper triangle of a matrix into a vector
+!****************************************************
     implicit none
     integer,intent(in)  :: n
     real(wp),intent(in) :: a(n,n)
@@ -147,7 +152,9 @@ subroutine dsqtoh(n,a,b)
   end subroutine dsqtoh
 
   subroutine dhtosq(n,a,b)
-    !converts upper triangle vector into a symmetric matrix
+!*********************************************************
+!* converts upper triangle vector into a symmetric matrix
+!*********************************************************
     implicit none
     integer,intent(in)  :: n
     real(wp),intent(out) :: a(n,n)
@@ -166,33 +173,47 @@ subroutine dhtosq(n,a,b)
     return
   end subroutine dhtosq
 
-  !Profjection of the translational and rotational contributions to the numerical Hessian plus the mass-weighting of the Hessian
+!=========================================================================================!
+
   subroutine prj_mw_hess(nat,at,nat3,xyz,hess)
+!***************************************************************
+!* Projection of the translational and rotational DOF out of 
+!* the numerical Hessian plus the mass-weighting of the Hessian
+!***************************************************************
     implicit none
 
     integer,intent(in) :: nat,nat3
-    integer :: at(nat),ich
+    integer :: at(nat)
     real(wp),intent(inout) :: hess(nat3,nat3)
-    real(wp) ::  hess_ut(nat3*(nat3+1)/2),pmode(nat3,1)
     real(wp) ::  xyz(3,nat)
+    !real(wp) ::  hess_ut(nat3*(nat3+1)/2),pmode(nat3,1)
+    real(wp),allocatable ::  hess_ut(:),pmode(:,:)
 
-    !Transforms matrix of the upper triangle vector
+    allocate(hess_ut(nat3*(nat3+1)/2), source=0.0_wp)
+    allocate(pmode(nat3,1), source=0.0_wp)
+
+    !> Transforms matrix of the upper triangle vector
     call dsqtoh(nat3,hess,hess_ut)
 
-    !Projection
+    !> Projection
     call trproj(nat,nat3,xyz,hess_ut,.false.,0,pmode,1)
 
-    !Transforms vector of the upper triangle into matrix
+    !> Transforms vector of the upper triangle into matrix
     call dhtosq(nat3,hess,hess_ut)
 
-    !Mass weighting
+    !> Mass weighting
     call mass_weight_hess(nat,at,nat3,hess)
 
+    deallocate(pmode,hess_ut)
   end subroutine prj_mw_hess
 
-  ! Prints the vibration spectrum of the a system. The intensity is only artficially included as 1000 for every vibration!!
-  subroutine print_vib_spectrum(nat,at,nat3,xyz,freq,dir,fname)
+!=========================================================================================!
 
+  subroutine print_vib_spectrum(nat,at,nat3,xyz,freq,dir,fname)
+!*********************************************************************
+!* Prints the frequencies in Turbomoles "vibspectrum" format 
+!* The intensity is only artficially set to 1000 for every vibration!!
+!**********************************************************************
     integer,intent(in) :: nat,nat3
     integer :: at(nat),i
     real(wp) ::  xyz(3,nat)
@@ -232,10 +253,13 @@ subroutine print_vib_spectrum(nat,at,nat3,xyz,freq,dir,fname)
 
   end subroutine print_vib_spectrum
 
-  !Prints the vibration spectrum of the a system as a g98.out.
-  !Routine is adapted from the xtb code.
-  subroutine print_g98_fake(nat,at,nat3,xyz,freq,hess,dir,fname)
+!=========================================================================================!
 
+  subroutine print_g98_fake(nat,at,nat3,xyz,freq,hess,dir,fname)
+!****************************************************************
+!* Prints the vibration spectrum of the a system as a g98.out.
+!* Routine is adapted from the xtb code.
+!****************************************************************
     integer,intent(in) :: nat,nat3
     integer :: at(nat)
     integer  :: gu,i,j,ka,kb,kc,la,lb,k
@@ -346,9 +370,13 @@ subroutine print_g98_fake(nat,at,nat3,xyz,freq,hess,dir,fname)
 
   end subroutine print_g98_fake
 
-  !Prints the numerical hessian
-  subroutine print_hessian(hess,nat3,dir,fname)
+!=========================================================================================!
 
+
+  subroutine print_hessian(hess,nat3,dir,fname)
+!*******************************
+!* Prints the numerical hessian
+!*******************************
     integer :: nat3,i,j,k
     real(wp) :: hess(nat3,nat3)
     character(len=*) :: fname
@@ -392,9 +420,13 @@ subroutine print_hessian(hess,nat3,dir,fname)
 
   end subroutine print_hessian
 
-  !Effective Hessian at an MECP is computed via Eq. 27 and Eq. 28 in https://doi.org/10.1002/qua.25124
-  subroutine effective_hessian(nat,nat3,grad1_i,grad2_i,hess1,hess2,heff)
+!=========================================================================================!
 
+  subroutine effective_hessian(nat,nat3,grad1_i,grad2_i,hess1,hess2,heff)
+!******************************************************************
+!* Effective Hessian at an MECP is computed via Eq. 27 and Eq. 28 
+!* in https://doi.org/10.1002/qua.25124
+!******************************************************************
     implicit none
     integer,intent(in) :: nat,nat3
     integer :: i,j,ii

src/algos/numhess.f90

@@ -116,8 +116,10 @@ subroutine crest_numhess(env,tim)
     allocate (hess(nat3,nat3,calc%ncalculations),source=0.0_wp)
     allocate (freq(nat3,n_freqs),source=0.0_wp)
 
-    !>-- Computes numerical Hessians and stores them individually for each level
+!*********************************************************************************
+!>--- Computes numerical Hessians and stores them individually for each level
     call numhess2(mol%nat,mol%at,mol%xyz,calc,hess,io)
+!*********************************************************************************
 
     write (stdout,*) 'done.'
     write (stdout,*)
@@ -179,6 +181,9 @@ subroutine crest_numhess(env,tim)
         !>-- Prints Hessian
         call print_hessian(hess(:,:,i),nat3,'','numhess'//trim(adjustl(atmp)))
 
+        !>--- Print dipole gradients (if they exist)
+        call calc%calcs(i)%dumpdipgrad('dipgrad'//trim(adjustl(atmp)))
+
         !>-- Projects and mass-weights the Hessian
         call prj_mw_hess(mol%nat,mol%at,nat3,mol%xyz,hess(:,:,i))
 

src/algos/parallel.f90

@@ -646,7 +646,7 @@ subroutine crest_search_multimd_init(env,mol,mddat,nsim)
     call defaultGF(env)
     write (stdout,*) 'list of applied metadynamics Vbias parameters:'
     do i = 1,env%nmetadyn
-      write (stdout,'(''$metadyn '',f10.5,f8.3,i5)') env%metadfac(i)/env%rednat,env%metadexp(i)
+      write (stdout,'(''$metadyn '',f10.5,f8.3,i5)') env%metadfac(i),env%metadexp(i)
     end do
     write (stdout,*)
 
@@ -695,7 +695,7 @@ subroutine crest_search_multimd_init2(env,mddats,nsim)
   allocate (mtds(nsim))
   do i = 1,nsim
     if (mddats(i)%simtype == type_mtd) then
-      mtds(i)%kpush = env%metadfac(i)/env%rednat
+      mtds(i)%kpush = env%metadfac(i)
       mtds(i)%alpha = env%metadexp(i)
       mtds(i)%cvdump_fs = float(env%mddump)
       mtds(i)%mtdtype = cv_rmsd

src/algos/singlepoint.f90

@@ -34,7 +34,7 @@ subroutine crest_singlepoint(env,tim)
   use crest_data
   use crest_calculator
   use strucrd
-  use gradreader_module, only: write_engrad
+  use gradreader_module,only:write_engrad
   implicit none
   type(systemdata),intent(inout) :: env
   type(timer),intent(inout)      :: tim
@@ -46,7 +46,7 @@ subroutine crest_singlepoint(env,tim)
   type(calcdata) :: calc
   real(wp) :: accuracy,etemp
 
-  real(wp) :: energy
+  real(wp) :: energy,dip
   real(wp),allocatable :: grad(:,:)
 
   character(len=*),parameter :: partial = '∂E/∂'
@@ -107,7 +107,7 @@ subroutine crest_singlepoint(env,tim)
         write (stdout,'(a12,a12,a10)') 'Atom A','Atom B','BO(A-B)'
         do i = 1,mol%nat
           do j = 1,i-1
-            if (calc%calcs(k)%wbo(i,j) > 0.0002_wp) then
+            if (calc%calcs(k)%wbo(i,j) > 0.01_wp) then
               write (stdout,*) i,j,calc%calcs(k)%wbo(i,j)
             end if
           end do
@@ -120,16 +120,33 @@ subroutine crest_singlepoint(env,tim)
     write (stdout,*)
   end if
 
-  if(all(calc%calcs(:)%rdgrad .eqv. .false.))then
+  if (any(calc%calcs(:)%rddip)) then
+    write (stdout,*)
+    write (stdout,*) 'Molecular dipole moments (a.u.):'
+    do k = 1,calc%ncalculations
+      if (calc%calcs(k)%rddip) then
+        dip = norm2(calc%calcs(k)%dipole)
+        write (stdout,'("> ",a,i0)') 'Calculation level ',k
+        write (stdout,'(a10,a10,a10,a12)') 'x','y','z','tot (Debye)'
+        write (stdout,'(4f10.3)') calc%calcs(k)%dipole(:),dip*autod
+      end if
+    end do
+    write (stdout,*)
+    write (stdout,'(a)') repeat('-',80)
+  else
+    write (stdout,*)
+  end if
+
+  if (all(calc%calcs(:)%rdgrad.eqv..false.)) then
     write (stdout,'(a)') '> No gradients calculated'
-  else 
-  write (stdout,'(a)') '> Final molecular gradient ( Eh/a0 ):'
-  write (stdout,'(13x,a,13x,a,13x,a)') partial//'x',partial//'y',partial//'z'
-  do i = 1,mol%nat
-    write (stdout,'(3f18.8)') grad(1:3,i)
-  end do
-  write (stdout,'(a,f18.8,a)') '> Gradient norm:',norm2(grad),' Eh/a0'
-  endif
+  else
+    write (stdout,'(a)') '> Final molecular gradient ( Eh/a0 ):'
+    write (stdout,'(13x,a,13x,a,13x,a)') partial//'x',partial//'y',partial//'z'
+    do i = 1,mol%nat
+      write (stdout,'(3f18.8)') grad(1:3,i)
+    end do
+    write (stdout,'(a,f18.8,a)') '> Gradient norm:',norm2(grad),' Eh/a0'
+  end if
 
   if (calc%ncalculations > 1) then
     write (stdout,*)
@@ -148,12 +165,12 @@ subroutine crest_singlepoint(env,tim)
   write (stdout,'(1x,a,f20.10,a)') 'GRADIENT NORM',norm2(grad),' Eh/a0'
   write (stdout,'(a)') repeat('=',40)
 
-  write(stdout,'(1x,a)') 'Writing crest.engrad ...' 
+  write (stdout,'(1x,a)') 'Writing crest.engrad ...'
   call write_engrad('crest.engrad',energy,grad)
 
-  if(env%testnumgrad)then
+  if (env%testnumgrad) then
     call numgrad(mol,calc,grad)
-  endif 
+  end if
 
   deallocate (grad)
 !========================================================================================!
@@ -173,7 +190,7 @@ subroutine crest_xtbsp(env,xtblevel,molin)
 !*  xtblevel - quick selection of calc. level
 !*  molin    - molecule data
 !********************************************************************
-  use crest_parameters 
+  use crest_parameters
   use crest_data
   use crest_calculator
   use strucrd
@@ -245,7 +262,7 @@ subroutine crest_ensemble_singlepoints(env,tim)
   use crest_calculator
   use strucrd
   use optimize_module
-  use utilities, only: dumpenergies
+  use utilities,only:dumpenergies
   implicit none
   type(systemdata),intent(inout) :: env
   type(timer),intent(inout)      :: tim
@@ -267,7 +284,7 @@ subroutine crest_ensemble_singlepoints(env,tim)
   real(wp) :: percent
   character(len=52) :: bar
 !========================================================================================!
- write (*,*)
+  write (*,*)
 !>--- check for the ensemble file
   inquire (file=env%ensemblename,exist=ex)
   if (ex) then
@@ -287,7 +304,7 @@ subroutine crest_ensemble_singlepoints(env,tim)
   call rdensemble(ensnam,nat,nall,at,xyz,eread)
 !>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<!
 !>--- Important: crest_oloop requires coordinates in Bohrs
-  xyz = xyz / bohr
+  xyz = xyz/bohr
 !>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<!
 
 !>--- set OMP parallelization
@@ -308,14 +325,14 @@ subroutine crest_ensemble_singlepoints(env,tim)
 
 !>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<!
 !>--- Important: ensemble file must be written in AA
-  xyz = xyz / angstrom
+  xyz = xyz/angstrom
 !>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<!
 !>--- write output ensemble
   call wrensemble(ensemblefile,nat,nall,at,xyz,eread)
-  write(stdout,'(/,a,a,a)') 'Ensemble with updated energies written to <',ensemblefile,'>'
+  write (stdout,'(/,a,a,a)') 'Ensemble with updated energies written to <',ensemblefile,'>'
 
   call dumpenergies('crest.energies',eread)
-  write(stdout,'(/,a,a,a)') 'List of energies written to <','crest.energies','>' 
+  write (stdout,'(/,a,a,a)') 'List of energies written to <','crest.energies','>'
 
   deallocate (eread,at,xyz)
 !========================================================================================!