fixes (start v3.0.1)

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)
 !========================================================================================!