SDF update (crest-lab#142)

* Strucreader changes from devel branch * Added SDF molfile conversion for protonation (etc.) tool ensembles * Further SDF handling cleanup * Some other minor changes and github workflow update
grimme-lab · Sep 27, 2022 · 19a8f00 · 19a8f00
1 parent 205b511
commit 19a8f00
Show file tree

Hide file tree

Showing 14 changed files with 1,730 additions and 1,598 deletions.
diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
@@ -3,7 +3,7 @@ on: [push, pull_request]
 
 env:
   BUILD_DIR: _build
-  XTB_VERSION: 6.4.0
+  XTB_VERSION: 6.5.0
 
 jobs:
   intel-meson-build:
@@ -49,7 +49,7 @@ jobs:
         printenv >> $GITHUB_ENV
 
     - name: Install meson
-      run: pip3 install meson==0.55.3 ninja
+      run: pip3 install meson==0.60.1 ninja
 
     - name: Configure meson build
       run: meson setup ${{ env.BUILD_DIR }} --prefix=/
@@ -73,14 +73,14 @@ jobs:
         owner: grimme-lab
         repo: xtb
         tag: v${{ env.XTB_VERSION }}
-        file: '/xtb-\d+\.tar\.xz/'
+        file: '/xtb-${{ env.XTB_VERSION }}-linux-x86_64\.tar\.xz/'
 
     - name: Add xtb to path
       run: |
         tar xvf xtb-*.tar.xz
         echo ${{ env.PWD }}/xtb-${{ env.XTB_VERSION }}/bin >> $GITHUB_PATH
 
-    - name: Run example 0
-      run: |
-        bash run.sh
-      working-directory: examples/expl-0
+#    - name: Run example 0
+#      run: |
+#        bash run.sh
+#      working-directory: examples/expl-0
diff --git a/src/Makefile b/src/Makefile
@@ -1,4 +1,4 @@
-PROG = ~/bin/crest
+PROG = ~/bin/crest/devel/crest
 
 OBJDIR = build
 #--------------------------------------------------------------------------

diff --git a/src/classes.f90 b/src/classes.f90
@@ -28,7 +28,6 @@ module crest_data
 
    public :: systemdata
    public :: timer
-   public :: sdfobj
    public :: protobj
    public :: constra
    public :: optlevflag
@@ -159,22 +158,6 @@ module crest_data
       logical,allocatable :: blacklist(:) !a blacklist of atoms to disallow deprotonation from
    end type protobj
 
-!-----------------------------------------------------------------------------------------------------
-!-----------------------------------------------------------------------------------------------------
-
-   type :: sdfobj
-      logical :: v3000 = .false. 
-      integer :: nat
-      integer :: nmisc
-      character(len=128) :: countsline
-      character(len=128),allocatable :: hblock(:)  !sdf header block (3 lines + counts line)
-      character(len=128),allocatable :: cblock(:)  !coordinate block (nat lines)
-      character(len=128),allocatable :: miscblock(:)  !misc block    (until-EOF lines)
-    contains
-      procedure :: deallocate => deallocate_sdf
-   end type sdfobj
-
-
 !-----------------------------------------------------------------------------------------------------
 !-----------------------------------------------------------------------------------------------------
 
@@ -368,9 +351,6 @@ module crest_data
     !--- saved constraints
       type(constra) :: cts
 
-    !--- SDF input format object
-      type(sdfobj) :: sdf
-
     !--- NCI mode data
       real(wp) :: potscal
 
@@ -480,6 +460,7 @@ module crest_data
       logical :: noconst=.false.   ! no constrain of solute during QCG Growth
       logical :: onlyZsort         ! do only the ZSORT routine ?
       logical :: optpurge = .false. !MDOPT purge application
+      logical :: outputsdf = .false. ! write output ensemble as sdf?
       logical :: pcaexclude = .false. ! exclude user set atoms from PCA?
       logical :: pclean            ! cleanup option for property mode
       logical :: performCross      ! perform the GC in V1/V2 ?
@@ -507,7 +488,7 @@ module crest_data
       logical :: scallen           ! scale the automatically determined MD length by some factor?
       logical :: scratch           ! use scratch directory
       logical :: setgcmax = .false.! adjust the maxmimum number of structures taken into account for GC?
-      logical :: sdfformat         ! was the SDF format used as input file?
+      logical :: sdfformat = .false. ! was the SDF format used as input file?
       logical :: slow              ! slowmode (counterpart to quick mode)
       logical :: solv_md = .false. !switches on QCG-ensemblerun instead of CFF
       logical :: staticmtd = .false. ! do a static MTD instead of normal MDs
@@ -654,13 +635,6 @@ subroutine stop_timer(self,n)
    self%t(n,3)=self%t(n,3) + (self%t(n,2)-self%t(n,1))
 end subroutine stop_timer
 !----------------------------------------------------------------------------------------------------
-subroutine deallocate_sdf(self)
-   implicit none
-   class(sdfobj) :: self
-   if (allocated( self%hblock)) deallocate( self%hblock )
-   if (allocated( self%cblock)) deallocate( self%cblock )
-   if (allocated( self%miscblock)) deallocate( self%miscblock )
-end subroutine deallocate_sdf
 !----------------------------------------------------------------------------------------------------
 subroutine allocate_file(self,n)
    implicit none

diff --git a/src/confparse.f90 b/src/confparse.f90
@@ -1292,6 +1292,10 @@ subroutine parseflags(env,arg,nra)
         if (io == 0 .and. (index(arg(i + 1),'-') .eq. 0)) then
           env%maxflip = nint(rdum)
         end if
+      case ('-osdf')
+        env%outputsdf = .true.
+        write(*,'(2x,a," :",1x,a)') trim(arg(i)), &
+        & "output ensemble requested in sdf format"
 !========================================================================================!
 !------ flags for parallelization / disk space
 !========================================================================================!
@@ -1827,7 +1831,8 @@ subroutine parseflags(env,arg,nra)
       case ('-pclean')                           !cleanup option for property mode, i.e., remove PROP/
         env%pclean = .true.
 !========================================================================================!
-      case ('-scratch')                          !use a scratch directory to perform the calculation in
+      case ('-scratch') 
+        !use a scratch directory to perform the calculation in
         env%scratch = .true.
         atmp = ''
         if (nra .ge. (i + 1)) atmp = adjustl(arg(i + 1))
@@ -2227,8 +2232,10 @@ subroutine inputcoords(env,arg)
   !>--- if the input was a SDF file, special handling
   env%sdfformat = .false.
   call checkcoordtype(inputfile,i)
-  if (i == 31 .or. i == 32) then
-    call inpsdf(env,inputfile)
+  if (any((/31,32/) == i)) then
+    !call inpsdf(env,inputfile)
+     env%sdfformat = .true.
+     env%outputsdf = .true.
   end if
 
   !>--- after this point there should always be an coord file present

diff --git a/src/confscript2_misc.f90 b/src/confscript2_misc.f90
@@ -35,8 +35,8 @@ subroutine xtbsp(env,xtblevel)
          !pipe=' > xtb.out 2>/dev/null'
          call remove('gfnff_topo')
          call remove('energy')
+         call remove('wbo')
          if(.not.env%chargesfile)call remove('charges')
-         call remove('xtbrestart')
 !---- (OPTIONAL) select xtb level and set flag
          if(present(xtblevel))then
              select case(xtblevel)
@@ -63,8 +63,6 @@ subroutine xtbsp(env,xtblevel)
          call copy('coord',fname)
          call clear_setblock(fname)
 !---- jobcall
-!         write(jobcall,'(a,1x,a,1x,a," --sp ",a)') &
-!         &     trim(env%ProgName),trim(fname),trim(xtbflag),trim(env%solv)
          jobcall = ""
          jobcall = trim(jobcall)//trim(env%ProgName)
          jobcall = trim(jobcall)//" "//trim(fname)//" --sp"
@@ -94,6 +92,7 @@ subroutine xtbsp2(fname,env)
          type(systemdata) :: env
          character(len=512) :: jobcall
          character(*),parameter :: pipe=' > xtbcalc.out 2>/dev/null'
+         character(len=4) :: chrgstr
          integer :: io
          call remove('gfnff_topo')
          call remove('energy')
@@ -104,9 +103,20 @@ subroutine xtbsp2(fname,env)
             call ompautoset(env%threads,7,env%omp,env%MAXRUN,1) !set the global OMP/MKL variables for the xtb jobs
          endif
 !---- jobcall
-         write(jobcall,'(a,1x,a,1x,a,'' --sp --wbo '',a,1x,a,a)') &
-         &     trim(env%ProgName),trim(fname),trim(env%gfnver),trim(env%solv),trim(pipe)
-         !call system(trim(jobcall))
+         jobcall = ""
+         jobcall = trim(jobcall)//trim(env%ProgName)
+         jobcall = trim(jobcall)//" "//trim(fname)//" --sp --wbo"
+         jobcall = trim(jobcall)//" "//trim(env%gfnver)
+         jobcall = trim(jobcall)//" "//trim(env%solv)
+         if( env%chrg /= 0 )then
+           write(chrgstr,'(i0)') env%chrg
+           jobcall = trim(jobcall)//" --chrg "//trim(chrgstr)
+         endif
+         if( env%uhf /= 0 )then
+           write(chrgstr,'(i0)') env%uhf
+           jobcall = trim(jobcall)//" --uhf "//trim(chrgstr)
+         endif
+         jobcall = trim(jobcall)//pipe
          call execute_command_line(trim(jobcall), exitstat=io)
 !---- cleanup
          call remove('xtbcalc.out')

diff --git a/src/crest_main.f90 b/src/crest_main.f90
@@ -155,29 +155,14 @@ program CREST
   !>--- protonation tool
        case( p_protonate )
         call protonate(env,tim)
-        if(env%relax)then
-          env%chrg=env%chrg + 1
-          env%nat = env%nat + 1
-          env%rednat = env%rednat +1
-          call relaxensemble('protonated.xyz',env,tim)
-        endif  
         call propquit(tim)
   !>--- deprotonation
        case( p_deprotonate )
         call deprotonate(env,tim)
-        if(env%relax)then
-           env%chrg=env%chrg-1
-           env%nat = env%nat - 1
-           env%rednat = env%rednat - 1
-           call relaxensemble('deprotonated.xyz',env,tim)
-        endif
         call propquit(tim)
   !>--- tautomerization
        case( p_tautomerize )
         call tautomerize(env,tim)
-        if(env%relax)then
-           call relaxensemble('tautomers.xyz',env,tim)
-        endif
         call propquit(tim)
   !>--- extended tautomerization
        case( p_tautomerize2 )
@@ -280,8 +265,13 @@ program CREST
            continue
       end select
 
-      if(env%sdfformat.and.(any((/1,2,22/)==env%crestver)))then    
-         call wrsdfens(env%sdf,conformerfile,conformerfilebase//'.sdf')
+      if(env%outputsdf .or. env%sdfformat) then
+        if(any((/crest_mfmdgc,crest_imtd,crest_imtd2/)==env%crestver))then    
+         call new_wrsdfens(env,conformerfile,conformerfilebase//'.sdf',.false.)
+        endif
+        if(any((/crest_screen,crest_mdopt/)==env%crestver))then
+         call new_wrsdfens(env,'crest_ensemble.xyz','crest_ensemble.sdf',.false.) 
+        endif
       endif
 
 !=========================================================================================!

diff --git a/src/deprotonate.f90 b/src/deprotonate.f90
@@ -161,6 +161,16 @@ subroutine deprotonate(env,tim)
      call sort_ens(deprot,'deprotonated.xyz',.true.)
      call tim%stop(2)
 
+!>--- (optional) post-processing
+     if(env%relax)then
+       env%rednat = env%rednat - 1
+       call relaxensemble('deprotonated.xyz',env,tim)
+     endif
+
+     if(env%outputsdf)then
+     call new_wrsdfens(env,'deprotonated.xyz','deprotonated.sdf',.true.)
+     endif
+
 !--- reset data for main dir
      env%chrg = refchrg
      if(env%chrg .eq. 0) then

diff --git a/src/printouts.f90 b/src/printouts.f90
@@ -22,7 +22,7 @@
 !===================================================================================================!
 subroutine confscript_head
       implicit none
-      character(len=40),parameter:: date='  Mo 20. Jun 12:32:33 BST 2022'
+      character(len=40),parameter:: date='Tue 27. Sep 01:07:11 BST 2022'
       character(len=10),parameter:: version='2.12  '
       logical :: niceprint
 
@@ -38,7 +38,7 @@ subroutine confscript_head
       write(*,'(7x,''|          Universitaet Bonn, MCTC           |'')')
       write(*,'(7x,''=============================================='')')
       write(*,'(7x,''Version '',a,'', '',a)')trim(version),trim(date)
-      write(*,'(2x,''Using the xTB program. Compatible with xTB version 6.4.0'')')
+      write(*,'(1x,''Using the xTB program. Compatible with xTB versions >6.4.0'')')
       write(*,*)
       write(*,'(3x,''Cite work conducted with this code as'')')
       write(*,'(/,3x,''• P.Pracht, F.Bohle, S.Grimme, PCCP, 2020, 22, 7169-7192.'')')

diff --git a/src/protonate.f90 b/src/protonate.f90
@@ -130,7 +130,7 @@ subroutine protonate(env,tim)
         call rename('OPTIM'//'/'//'opt.xyz',trim(filename))
         call rmrf('OPTIM')
         !write(*,*) trim(filename)
-        call copy(trim(filename),'ensemble-test.xyz')
+        !call copy(trim(filename),'ensemble-test.xyz')
         if(prot%ABcorrection)then
             call prot_correction(env,trim(filename))
         endif
@@ -173,6 +173,19 @@ subroutine protonate(env,tim)
      call sort_ens(prot,'protonated.xyz',.true.)
      call tim%stop(2)
 
+
+
+!>--- (optional) post-processing
+     if(env%relax)then
+       env%rednat = env%rednat +1
+       call relaxensemble('protonated.xyz',env,tim)
+     endif
+
+     if(env%outputsdf)then
+     call new_wrsdfens(env,'protonated.xyz','protonated.sdf',.true.)
+     endif  
+
+
 !--- reset data for main dir
      env%chrg = refchrg
      if(env%chrg .eq. 0) then
@@ -183,7 +196,6 @@ subroutine protonate(env,tim)
        close(ich)
      endif
      env%nat=natp - 1 !reset nat
-     !call chdir(thispath)
 end subroutine protonate
 
 !--------------------------------------------------------------------------------------------
@@ -195,11 +207,10 @@ subroutine xtblmo(env)
          use crest_data
          implicit none
          type(systemdata) :: env
-         character(len=80) :: fname,pipe
+         character(len=80) :: fname
          character(len=512) :: jobcall
          integer :: io
-!---- some options
-         pipe=' > xtb.out 2>/dev/null'
+         character(len=*),parameter :: pipe = ' > xtb.out 2>/dev/null'
 
 !---- setting threads
          if(env%autothreads)then
@@ -212,8 +223,9 @@ subroutine xtblmo(env)
 !---- jobcall
          write(*,*)
          write(*,'('' LMO calculation ... '')',advance='no')
-         write(jobcall,'(a,1x,a,1x,a,'' --sp --lmo '',a,1x,a,a)') &
-         &     trim(env%ProgName),trim(fname),trim(env%gfnver),trim(env%solv),trim(pipe)
+         write(jobcall,'(a,1x,a,1x,a,'' --sp --lmo'',1x,a)') &
+         &     trim(env%ProgName),trim(fname),trim(env%gfnver),trim(env%solv)
+         jobcall = trim(jobcall)//trim(pipe)
          !call system(trim(jobcall))
          call execute_command_line(trim(jobcall), exitstat=io)
          write(*,'(''done.'')')

diff --git a/src/qcg/solvtool_misc.f90 b/src/qcg/solvtool_misc.f90
@@ -738,7 +738,7 @@ subroutine wr_cluster_cut(fname_cluster,n1,n2,iter,fname_solu_cut,fname_solv_cut
   character(len=*),intent(in) :: fname_cluster, fname_solu_cut,fname_solv_cut
   character (len=256)         :: atmp
   character (len=2)           :: a2   
-  integer                     :: ich,i,k,stat,io
+  integer                     :: ich,i,k,stat,io,io2
 
 
   ich=142
@@ -749,7 +749,7 @@ subroutine wr_cluster_cut(fname_cluster,n1,n2,iter,fname_solu_cut,fname_solv_cut
      read(ich,'(a)',iostat=io) atmp
      if(io < 0) exit
        atmp = adjustl(atmp) 
-       call coordline(atmp,a2,xyz1(1:3,k))
+       call coordline(atmp,a2,xyz1(1:3,k),io2)
        at1(k) = e2i(a2)
      k=k+1
   end do
@@ -758,7 +758,7 @@ subroutine wr_cluster_cut(fname_cluster,n1,n2,iter,fname_solu_cut,fname_solv_cut
      read(ich,'(a)',iostat=io) atmp
      if(io < 0) exit
        atmp = adjustl(atmp) 
-       call coordline(atmp,a2,xyz2(1:3,k))
+       call coordline(atmp,a2,xyz2(1:3,k),io2)
        at2(k) = e2i(a2)
      k=k+1
   end do