Fix QCG constraining & Cleanup

Signed-off-by: cplett <82893466+cplett@users.noreply.github.com>

src/qcg/solvtool.f90

@@ -170,7 +170,7 @@ subroutine qcg_setup(env, solu, solv)
    type(systemdata):: env
    type(zmolecule) :: solv, solu
 
-   integer :: io, f, r, v, ich
+   integer :: io, f, r
    integer :: num_O, num_H, i
    character(len=*), parameter :: outfmt = '(1x,1x,a,1x,f14.7,a,1x)'
    logical :: e_there, tmp, used_tmp
@@ -179,41 +179,26 @@ subroutine qcg_setup(env, solu, solv)
    character(len=80)  :: atmp
    character(len=20)  :: gfnver_tmp
 
-
    call getcwd(thispath)
 
+   ! Remove scratch dir, if present
    inquire (file='./qcg_tmp/solute_properties/solute', exist=tmp)
    if (tmp) call rmrf('qcg_tmp') !User given scratch dir will be removed anyway after run
 
-   if (env%scratchdir .eq. '') then !check if scratch was not set
+   !Make scratch directories
+   if (env%scratchdir .eq. '') then
       env%scratchdir = 'qcg_tmp'
       io = makedir('qcg_tmp')
    end if
-
-   call copysub('solute', trim(env%scratchdir))
-   call copysub('solvent', env%scratchdir)
-   call env%wrtCHRG(env%scratchdir) !Write .CHRG and .UHF with 3 lines for xtb docking
+   if(env%fixfile /= 'none selected') then
+      call copysub(env%fixfile, env%scratchdir)
+   end if
    call chdir(env%scratchdir)
-
    f = makedir('solute_properties')
+   if(env%fixfile /= 'none selected') then
+      call copysub(env%fixfile, env%scratchdir)
+   end if
    r = makedir('solvent_properties')
-   v = makedir('tmp_grow')
-
-   call chdir('tmp_grow')
-   call getcwd(tmp_grow)
-   call chdir(thispath)
-   call chdir(env%scratchdir)
-
-   call copysub('solute', 'solute_properties')
-   call copysub('.CHRG', 'solute_properties')
-   call copysub('.UHF', 'solute_properties')
-   call copysub('solvent', 'solvent_properties')
-
-   call copysub('.CHRG', tmp_grow)
-   call copysub('.UHF', tmp_grow)
-
-   call remove('solute')
-   call remove('solvent')
 
    if (.not. env%nopreopt) then
       write (*, *)
@@ -227,9 +212,7 @@ subroutine qcg_setup(env, solu, solv)
 
 !---- Properties solute
    call chdir('solute_properties')
-   env%ref%nat = solu%nat
-   env%ref%at = solu%at
-   env%ref%xyz = solu%xyz
+   call env%wrtCHRG('') !Write three lines in QCG mode, but xtb anyway only reads first one
 
 !---- Geometry preoptimization solute
    if (env%final_gfn2_opt) then !If GFN2 final opt, solute also GFN2 optimized
@@ -238,31 +221,7 @@ subroutine qcg_setup(env, solu, solv)
    end if
 
    if ((.not. env%nopreopt) .and. (solu%nat /= 1)) then
-      call wrc0('coord', solu%nat, solu%at, solu%xyz) !write coord for xtbopt routine
-      if (env%cts%used) then
-         call wrc0('coord.ref', solu%nat, solu%at, solu%xyz) !write coord for xtbopt routine
-      end if
-!---- coord setup section
-      open (newunit=ich, file='coord')
-      do
-         read (ich, '(a)', iostat=io) atmp
-         if (io < 0) exit
-         if (index(atmp, '$coord') .ne. 0) cycle
-         if (index(atmp(1:1), '$') .ne. 0) then
-            !write(ich,'(a)')'$end'
-            exit
-         end if
-      end do
-      !add constraints (only if given, else the routine returns)
-      call write_cts(ich, env%cts)
-      call write_cts_biasext(ich, env%cts)
-      write (ich, '(a)') '$end'
-      close (ich)
-
-      !call xtbopt(env)
-      call trialOPT(env)
-      call rdcoord('coord', solu%nat, solu%at, solu%xyz)
-      call remove('coord')
+      call xtb_opt_qcg(env, solu, .true.)
    end if
 
 !--- Axistrf
@@ -272,9 +231,9 @@ subroutine qcg_setup(env, solu, solv)
 !---- LMO/SP-Computation solute
    if (env%use_xtbiff) then
       write (*, *) 'Generating LMOs for solute'
-      call xtb_lmo(env, 'solute')!,solu%chrg)
+      call xtb_lmo(env, 'solute')
    else
-      call xtbsp3(env, 'solute')
+      call xtb_sp_qcg(env, 'solute')
    end if
 
    if (env%final_gfn2_opt) then !If GFN2 final opt, solute also GFN2 optimized
@@ -290,9 +249,6 @@ subroutine qcg_setup(env, solu, solv)
 
    if (env%use_xtbiff) then
       call rename('xtblmoinfo', 'solute.lmo')
-      call copysub('solute.lmo', tmp_grow)
-   else
-      call copysub('solute', tmp_grow)
    end if
 
    call chdir(thispath)
@@ -304,17 +260,11 @@ subroutine qcg_setup(env, solu, solv)
 !---- Properties solvent
    call chdir(env%scratchdir)
    call chdir('solvent_properties')
-   env%ref%nat = solv%nat
-   env%ref%at = solv%at
-   env%ref%xyz = solv%xyz
+   !No charges for solvent written. This is currently not possible
 
 !---- Geometry preoptimization solvent
    if ((.not. env%nopreopt) .and. (solv%nat /= 1)) then
-      call wrc0('coord', solv%nat, solv%at, solv%xyz) !write coord for xtbopt routine
-      !call xtbopt(env)
-      call trialOPT(env)
-      call rdcoord('coord', solv%nat, solv%at, solv%xyz)
-      call remove('coord')
+      call xtb_opt_qcg(env, solv, .false.)
    end if
    call wrc0('solvent', solv%nat, solv%at, solv%xyz)
 
@@ -323,7 +273,7 @@ subroutine qcg_setup(env, solu, solv)
       write (*, *) 'Generating LMOs for solvent'
       call xtb_lmo(env, 'solvent')!,solv%chrg)
    else
-      call xtbsp3(env, 'solvent')
+     call xtb_sp_qcg(env, 'solvent')
    end if
 
    call grepval('xtb.out', '| TOTAL ENERGY', e_there, solv%energy)
@@ -335,9 +285,6 @@ subroutine qcg_setup(env, solu, solv)
 
    if (env%use_xtbiff) then
       call rename('xtblmoinfo', 'solvent.lmo')
-      call copysub('solvent.lmo', tmp_grow)
-   else
-      call copysub('solvent', tmp_grow)
    end if
 
    call chdir(thispath)
@@ -379,21 +326,14 @@ subroutine read_qcg_input(env, solu, solv)
    logical                        :: pr
    real(wp), parameter            :: amutokg = 1.66053886E-27
    real(wp), parameter            :: third = 1.0d0/3.0d0
-   integer                        :: i, ich
+   integer                        :: i
    real(wp)                       :: r_solu, r_solv
 
    pr = .true.
 
+!--- Read in solu and solv coordinates and make solute and solvent file in WD
    call inputcoords_qcg(env, solu, solv)
 
-!--- Read in nat, at, xyz
-!   call simpletopo_file('solute', solu, .false., .false., '')
-!   allocate (solu%xyz(3, solu%nat))
-!   call rdcoord('solute', solu%nat, solu%at, solu%xyz)
-!   call simpletopo_file('solvent', solv, .false., .false., '')
-!   allocate (solv%xyz(3, solv%nat))
-!   call rdcoord('solvent', solv%nat, solv%at, solv%xyz)
-
 !--- CMA-Trafo
    call cma_shifting(solu, solv)
 
@@ -520,7 +460,7 @@ subroutine qcg_grow(env, solu, solv, clus, tim)
 
    integer                    :: minE_pos, m
    integer                    :: iter = 1
-   integer                    :: i, j, io
+   integer                    :: i, j, io, v
    integer                    :: max_cycle
    logical                    :: e_there, high_e, success, neg_E
    real(wp)                   :: etmp(500)
@@ -604,10 +544,20 @@ end subroutine both_ellipsout
    call get_ellipsoid(env, solu, solv, clus, .true.)
    call pr_grow_energy()
 
-   if (env%cts%used) call copysub(env%solu_file, env%scratchdir)
+   !Set up directory
    call chdir(env%scratchdir)
-   if (env%cts%used) call copysub(env%solu_file, 'tmp_grow')
+   v = makedir('tmp_grow')
+   if(env%fixfile /= 'none selected') then
+      call copysub(env%fixfile, 'tmp_grow')
+   end if
+   if (env%use_xtbiff) then
+     call copy('solute_properties/solute.lmo', 'tmp_grow/solute.lmo')
+     call copy('solvent_properties/solvent.lmo', 'tmp_grow/solvent.lmo')
+   end if
+
    call chdir('tmp_grow')
+   call wrc0('solute', solu%nat, solu%at, solu%xyz)
+   call wrc0('solvent', solv%nat, solv%at, solv%xyz)
 
    call ellipsout('solute_cavity.coord', clus%nat, clus%at, clus%xyz, solu%ell_abc)
    solv%ell_abc = clus%ell_abc
@@ -657,7 +607,7 @@ end subroutine both_ellipsout
                call check_dock(neg_E)
             end if
 
-!!! If Interaction Energy is not negativ and existent, wall pot. too small and increase
+!-- If Interaction Energy is not negativ and existent, wall pot. too small and increase
             if (neg_E) then
                success = .true.
             else
@@ -1676,7 +1626,7 @@ end subroutine aver
       call chdir(to)
       call wrc0('cluster.coord', clus%nat, clus%at, clus%xyz)
       call wr_cluster_cut('cluster.coord', solu%nat, solv%nat, env%nsolv, 'solute_cut.coord', 'solvent_shell.coord')
-      call xtbsp3(env, 'solvent_shell.coord')
+      call xtb_sp_qcg(env, 'solvent_shell.coord')
       call grepval('xtb.out', '| TOTAL ENERGY', ex, e_empty(i))
       call copy('solvent_shell.coord', 'solvent_cluster.coord')
       call copy('solvent_cluster.coord', 'filled_cluster.coord')
@@ -2807,13 +2757,13 @@ subroutine get_interaction_E(env, solu, solv, clus, iter, E_inter)
    call wr_cluster_cut('cluster.coord', solu%nat, solv%nat, iter, 'solute_cut.coord', 'solvent_cut.coord')
 
 !--- Perform single point calculations and recieve energies
-   call xtbsp3(env, 'solute_cut.coord')
+   call xtb_sp_qcg(env, 'solute_cut.coord')
    call grepval('xtb.out', '| TOTAL ENERGY', e_there, e_solute)
    if (.not. e_there) write (*, *) 'Solute energy not found'
-   call xtbsp3(env, 'solvent_cut.coord')
+   call xtb_sp_qcg(env, 'solvent_cut.coord')
    call grepval('xtb.out', '| TOTAL ENERGY', e_there, e_solvent)
    if (.not. e_there) write (*, *) 'Solvent energy not found'
-   call xtbsp3(env, 'cluster.coord')
+   call xtb_sp_qcg(env, 'cluster.coord')
    call grepval('xtb.out', '| TOTAL ENERGY', e_there, e_cluster)
    if (.not. e_there) write (*, *) 'Cluster energy not found'
 

src/qcg/solvtool_misc.f90

@@ -20,7 +20,7 @@
 !--------------------------------------------------------------------------------------------
 ! A quick single point xtb calculation without wbo
 !--------------------------------------------------------------------------------------------
-subroutine xtbsp3(env, fname)
+subroutine xtb_sp_qcg(env, fname)
    use iso_fortran_env, only: wp => real64
    use iomod
    use crest_data
@@ -49,7 +49,57 @@ subroutine xtbsp3(env, fname)
    call remove('xtbrestart')
    call remove('xtbtopo.mol')
    call remove('gfnff_topo')
-end subroutine xtbsp3
+end subroutine xtb_sp_qcg
+
+!--------------------------------------------------------------------------------------------
+! A quick single xtb optimization gets zmol and overwrites it with optimized stuff
+!--------------------------------------------------------------------------------------------
+subroutine xtb_opt_qcg(env, zmol, constrain)
+   use iso_fortran_env, only: wp => real64
+   use iomod
+   use crest_data
+   use zdata
+   use strucrd
+
+   implicit none
+   type(systemdata), intent(in) :: env
+   type(zmolecule), intent(inout) :: zmol
+
+   character(:), allocatable :: fname
+   character(len=512) :: jobcall
+   logical :: constrain
+   logical :: const
+   character(*), parameter :: pipe = ' > xtb_opt.out 2> /dev/null'
+   integer :: io
+
+!--- Write coordinated
+   fname = 'coord'
+   call wrc0(fname, zmol%nat, zmol%at, zmol%xyz) !write coord for xtbopt routine
+
+!---- setting threads
+   if (env%autothreads) then
+      call ompautoset(env%threads, 7, env%omp, env%MAXRUN, 1) !set the global OMP/MKL variables for the xtb jobs
+   end if
+!---- jobcall & Handling constraints
+   if(constrain .AND. env%cts%used) then
+     call write_constraint(env, fname, 'xcontrol')
+     call wrc0('coord.ref', zmol%nat, zmol%at, zmol%xyz) !write coord for xtbopt routine
+     write (jobcall, '(a,1x,a,1x,a,'' --opt --input xcontrol '',a,1x,a)') &
+     &     trim(env%ProgName), trim(fname), trim(env%gfnver), trim(env%solv), trim(pipe)
+   else
+     write (jobcall, '(a,1x,a,1x,a,'' --opt '',a,1x,a)') &
+     &     trim(env%ProgName), trim(fname), trim(env%gfnver), trim(env%solv), trim(pipe)
+   end if
+
+   call command(trim(jobcall), io)
+!---- cleanup
+   call rdcoord('xtbopt.coord', zmol%nat, zmol%at, zmol%xyz)
+   call remove('energy')
+   call remove('charges')
+   call remove('xtbrestart')
+   call remove('xtbtopo.mol')
+   call remove('gfnff_topo')
+end subroutine xtb_opt_qcg
 
 !___________________________________________________________________________________
 !
@@ -143,6 +193,7 @@ subroutine xtb_dock(env, fnameA, fnameB, solu, clus)
    use iomod
    use crest_data
    use zdata
+   use strucrd
 
    implicit none
 
@@ -955,6 +1006,31 @@ subroutine check_iff(neg_E)
 
 end subroutine check_iff
 
+!----------------------------------------------------------------------------
+! write constraints in a file used as xtb input
+
+subroutine write_constraint(env,coord_name,fname)
+  use iso_fortran_env, only : wp => real64
+  use crest_data
+  use iomod
+
+  implicit none
+
+  type(systemdata)     :: env
+  character(len=*),intent(in)    :: fname, coord_name
+
+  call copysub(coord_name, 'coord.ref')
+  open(unit=31,file=fname)
+  call write_cts(31,env%cts)
+  call write_cts_biasext(31,env%cts)
+  if(env%cts%used) then !Only, if user set constrians is an $end written
+     write(31,'(a)') '$end'
+  end if
+
+  close(31)
+
+end subroutine write_constraint
+
 !----------------------------------------------------------------------------
 ! write a wall potential in a file used as xtb input