diff --git a/src/classes.f90 b/src/classes.f90 index 0466bb79..ca310b60 100644 --- a/src/classes.f90 +++ b/src/classes.f90 @@ -349,7 +349,7 @@ module crest_data character(len=512) :: homedir !> original directory from which calculation was started character(len=512) :: scratchdir !> path to the scratch directory character(len=:),allocatable :: cmd - character(len=:),allocatable :: inputcoords,inputcoords_solv,inputcoords_solu + character(len=:),allocatable :: inputcoords character(len=:),allocatable :: wbofile character(len=:),allocatable :: atlist character(len=:),allocatable :: chargesfilename @@ -385,8 +385,6 @@ module crest_data !>--- property data objects type(protobj) :: ptb - type(protobj) :: ptb_solvent - type(protobj) :: ptb_solute !>--- saved constraints type(constra) :: cts @@ -408,8 +406,6 @@ module crest_data !>--- reference structure data (the input structure) type(refdata) :: ref - type(refdata) :: qcg_solvent - type(refdata) :: qcg_solute !>--- QCG data integer :: qcg_runtype = 0 !> Default is grow, 1= ensemble & opt, 2= e_solv, 3= g_solv diff --git a/src/confparse.f90 b/src/confparse.f90 index 741e6fa0..bd6aa967 100644 --- a/src/confparse.f90 +++ b/src/confparse.f90 @@ -491,6 +491,12 @@ subroutine parseflags(env,arg,nra) case ('-solvtool','-qcg') ! env%mddumpxyz = 1000 + ! Set solute file if present + if(i == 2) env%solu_file = trim(arg(i-1)) + ! Set solvent file if prensent + ! If it is another argument, it doesent matter as solvent file is checke in solvtool + if (nra >= i+1) env%solv_file = trim(arg(i+1)) + ! Set QCG defaults env%preopt = .false. env%crestver = crest_solv env%QCG = .true. @@ -501,9 +507,6 @@ subroutine parseflags(env,arg,nra) env%ewin = 3.0d0 env%doOHflip = .false. !Switch off OH-flip if (env%iterativeV2) env%iterativeV2 = .false. - if (nra >= i+1) then - env%solv_file = trim(arg(i+1)) - end if exit case ('-compress') @@ -2369,8 +2372,7 @@ subroutine inputcoords(env,arg) !>--- Redirect for QCG input reading if (env%QCG) then - arg2 = env%solv_file - call inputcoords_qcg(env,arg,arg2) + !Input coordinates are processed in solvtool.f90 file during solvtool subroutine return end if !>--- @@ -2444,162 +2446,3 @@ subroutine inputcoords(env,arg) return end subroutine inputcoords - -!========================================================================================! -!> Convert given QCG coordinate files into (TM format) -!> If the input is in SDF format, document the info to convert the -!> final ensemble back into this format -!========================================================================================! -subroutine inputcoords_qcg(env,arg1,arg2) - use iso_fortran_env,only:wp => real64 - use crest_data - use strucrd - use zdata - use iomod - implicit none - - type(systemdata) :: env - character(len=*) :: arg1,arg2 - - logical :: ex11,ex12,ex21,ex22,solu,solv - character(len=:),allocatable :: inputfile - type(coord) :: mol - type(zmolecule) :: zmol,zmol1 - integer :: i - -!--------------------Checking for input-------------! - - inquire (file='solute',exist=solu) - inquire (file=arg1,exist=ex11) - inquire (file='coord',exist=ex12) - - inquire (file='solvent',exist=solv) - inquire (file=arg2,exist=ex21) - if (len_trim(arg2) .eq. 0) then !Check if the second argument is just empty - ex21 = .false. - end if - inquire (file='coord',exist=ex22) - - if (.not.ex11.and..not.ex12.and..not.solu) then - error stop 'No (valid) solute file! exit.' - else if (.not.ex21.and..not.ex22.and..not.solv) then - error stop 'No (valid) solvent file! exit.' - end if - -!---------------Handling solute---------------------! - - if (ex11.and.arg1(1:1) .ne. '-') then - call mol%open(arg1) - call mol%write('solute') - call mol%write('solute.xyz') - call mol%deallocate() - inputfile = trim(arg1) - write (*,*) 'Solute-file: ',arg1 - env%solu_file = arg1 - else if (solu.and..not.ex11) then - call copy('solute','solute.old') - inputfile = 'solute' - write (*,'(/,1x,a)') 'Solute-file: solute' - env%solu_file = 'solute' - else if (ex12.and..not.ex11) then !-- save coord as reference - call copy('coord','coord.old') - call copy('coord','solute') - write (*,'(/,1x,a)') 'Solute-file: coord' - inputfile = 'coord' - env%solu_file = 'coord' - else - write (*,*) 'An error occured processing the solute-input file' - end if - - !--- 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) - write (*,*) 'QCG currently does not support SDF-files.' - end if - -!---------------Handling solvent---------------------! - - if (ex21.and.arg2(1:1) .ne. '-'.and.arg2(1:1) .ne. ' ') then - call mol%open(arg2) - call mol%write('solvent') - call mol%write('solvent.xyz') - call mol%deallocate() - inputfile = trim(arg2) - write (*,*) 'Solvent-file: ',arg2 - env%solv_file = arg2 - else if (solv.and..not.ex21) then - call copy('solvent','solvent.old') - inputfile = 'solvent' - write (*,'(/,1x,a)') 'Solvent-file: solvent' - env%solv_file = 'solvent' - else if (ex22.and..not.ex21) then !-- save coord as reference - call copy('coord','coord.old') - call copy('coord','solvent') - inputfile = 'coord' - write (*,'(/,1x,a)') 'Solvent-file: coord' - env%solv_file = 'coord' - else - write (*,*) 'An error occured during solvent-file processing' - end if - - !--- 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) - write (*,*) 'QCG currently does not support SDF-files.' - end if - -!-------------Checking both files and saving them---------------------------! - - !--- after this point there should always be a solvent and solute file present - if (.not.allocated(env%inputcoords_solu)) env%inputcoords_solu = 'solute' - if (.not.allocated(env%inputcoords_solv)) env%inputcoords_solv = 'solvent' - - call mol%open('solute') - - env%nat = mol%nat - !--- solute geo - env%qcg_solute%nat = mol%nat - env%qcg_solute%at = mol%at - env%qcg_solute%xyz = mol%xyz - if (any((/crest_mfmdgc,crest_imtd,crest_imtd2/) == env%crestver)) then - if (.not.env%autozsort) then - env%qcg_solute%ntopo = mol%nat*(mol%nat+1)/2 - allocate (env%qcg_solute%topo(env%qcg_solute%ntopo)) - call quicktopo(mol%nat,mol%at,mol%xyz,env%qcg_solute%ntopo,env%qcg_solute%topo) - end if - end if - call mol%deallocate() - - call mol%open('solvent') - - env%nat = mol%nat - env%rednat = env%nat !get the number of atoms and the reduced number of atoms if some of them are excluded from the RMSD calc in V2 - !--- solvent geo - env%qcg_solvent%nat = mol%nat - env%qcg_solvent%at = mol%at - env%qcg_solvent%xyz = mol%xyz - if (any((/crest_mfmdgc,crest_imtd,crest_imtd2/) == env%crestver)) then - if (.not.env%autozsort) then - env%qcg_solvent%ntopo = mol%nat*(mol%nat+1)/2 - allocate (env%qcg_solvent%topo(env%qcg_solvent%ntopo)) - call quicktopo(mol%nat,mol%at,mol%xyz,env%qcg_solvent%ntopo,env%qcg_solvent%topo) - end if - end if - call mol%deallocate() - - !--- for protonation/deprotonation applications get ref. number of fragments - !--- also get some other structure based info - call simpletopo_file('solute',zmol,.false.,.false.,'') - env%ptb_solute%nfrag = zmol%nfrag - call zmol%deallocate() - - call simpletopo_file('solvent',zmol1,.false.,.false.,'') - env%ptb_solvent%nfrag = zmol1%nfrag - call zmol1%deallocate() - - return -end subroutine inputcoords_qcg diff --git a/src/crest_main.f90 b/src/crest_main.f90 index 1e3b0912..59c99abc 100644 --- a/src/crest_main.f90 +++ b/src/crest_main.f90 @@ -88,15 +88,7 @@ program CREST !> SOME I/O STUFF !=========================================================================================! !>--- check for the coord file in the working directory - if (env%crestver .eq. crest_solv) then - inquire (file='solute',exist=ex1) - inquire (file='solvent',exist=ex2) - if (.not.ex1) then - error stop 'No solute file found. Exit.' - else if (.not.ex2) then - error stop 'No solvent file found. Exit.' - end if - else + if (env%crestver /= crest_solv) then inquire (file='coord',exist=ex) if (.not.ex) then error stop 'No coord file found. Exit.' diff --git a/src/qcg/solvtool.f90 b/src/qcg/solvtool.f90 index b828347a..13234c4f 100644 --- a/src/qcg/solvtool.f90 +++ b/src/qcg/solvtool.f90 @@ -1,7 +1,7 @@ !================================================================================! ! This file is part of crest. ! -! Copyright (C) 2021 Sebastian Spicher, Christoph Plett, Philipp Pracht +! Copyright (C) 2023 Christoph Plett, Sebastian Spicher, Philipp Pracht ! ! crest is free software: you can redistribute it and/or modify it under ! the terms of the GNU Lesser General Public License as published by @@ -55,11 +55,12 @@ subroutine crest_solvtool(env, tim) !> Check, if xtb is present call checkprog_silent(env%ProgName,.true.,iostat=io) - if(io /= 0 ) error stop + if(io /= 0 ) error stop 'No xtb found' -!> Check, if xtbiff is present +!> Check, if xtbiff is present (if it is required) if (env%use_xtbiff) then - call check_prog_path_iff(env) + call checkprog_silent(env%ProgIFF,.true.,iostat=io) + if(io /= 0 ) error stop 'No xtbiff found' else write (*, *) write (*, *) ' The use of the aISS algorithm is requested (recommend).' @@ -383,13 +384,15 @@ subroutine read_qcg_input(env, solu, solv) pr = .true. + 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) +! 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) @@ -428,7 +431,6 @@ subroutine read_qcg_input(env, solu, solv) call read_directed_input(env) end if - end subroutine read_qcg_input !> Read input for directed docking @@ -3331,32 +3333,6 @@ subroutine qcg_cleanup(env) end subroutine qcg_cleanup -subroutine check_prog_path_iff(env) - use crest_data - use iomod, only: command,checkprog_silent - implicit none - type(systemdata):: env ! MAIN STORAGE OS SYSTEM DATA - character(len=512) :: prog - character(len=256) :: str - character(len=256) :: path - integer :: ios, io - - prog = env%ProgIFF - !write (str, '("which ",a," > ",a,"_path 2>/dev/null")') trim(prog), trim(prog) - !call command(trim(str), io) - !write (str, '(a,"_path")') trim(prog) - !str = trim(str) - !open (unit=27, file=str, iostat=ios) - !read (27, '(a)', iostat=ios) path - call checkprog_silent( prog, iostat=ios) - if (ios .ne. 0) then - write (0, *) 'No xtb-IFF found. This is currently required for ', & - & 'QCG and available at https:/github.com/grimme-lab/xtbiff/' - error stop - end if - -end subroutine check_prog_path_iff - subroutine write_reference(env, solu, clus) use iso_fortran_env, wp => real64 use crest_data @@ -3377,3 +3353,80 @@ subroutine write_reference(env, solu, clus) call wrc0(env%fixfile, ref_clus%nat, ref_clus%at, ref_clus%xyz) end subroutine write_reference + + +!========================================================================================! +!> Convert given QCG coordinate files into (TM format) +!> Write "solute" and "solvent" coordinate files +!========================================================================================! +subroutine inputcoords_qcg(env, solute, solvent) + use iso_fortran_env,only:wp => real64 + use crest_data + use strucrd + use zdata + use iomod + implicit none + + type(systemdata), intent(inout) :: env + type(zmolecule), intent(out) :: solute, solvent + + logical :: ex11,ex21,solu,solv + type(coord) :: mol + type(zmolecule) :: zmol,zmol1 + integer :: i + +!--------------------Checking for input-------------! + + !Solute + inquire (file=env%solu_file,exist=ex11) + inquire (file='solute',exist=solu) + if (solu) call copy('solute','solute.old') !Backup solute file + if ((.not. ex11) .and. (.not. solu)) then + error stop 'No (valid) solute file! exit.' + else if ((.not. ex11) .and. (solu)) then + env%solu_file = 'solute' + end if + + !Solvent + inquire (file=env%solv_file,exist=ex21) + inquire (file='solvent',exist=solv) + if (solu) call copy('solvent','solvent.old') !Backup solvent file + if ((.not. ex21) .and. (.not. solv)) then + error stop 'No (valid) solvent file! exit.' + else if ((.not. ex11) .and. (solu)) then + env%solu_file = 'solvent' + end if + +!---------------Handling solute---------------------! + call mol%open(env%solu_file) + call mol%write('solute') + solute%nat = mol%nat + solute%at = mol%at + solute%xyz = mol%xyz + call mol%deallocate() + + !--- if the input was a SDF file, special handling + env%sdfformat = .false. + call checkcoordtype(env%solu_file,i) + if (i == 31.or.i == 32) then + !Add sdf stuff here, if somebody needs it + end if + +!---------------Handling solvent---------------------! + + call mol%open(env%solv_file) + call mol%write('solvent') + solvent%nat = mol%nat + solvent%at = mol%at + solvent%xyz = mol%xyz + call mol%deallocate() + + !--- if the input was a SDF file, special handling + env%sdfformat = .false. + call checkcoordtype(env%solv_file,i) + if (i == 31.or.i == 32) then + !Add sdf stuff here, if somebody needs it + end if + + return +end subroutine inputcoords_qcg diff --git a/src/qcg/solvtool_misc.f90 b/src/qcg/solvtool_misc.f90 index f3cd2a92..8739cb97 100644 --- a/src/qcg/solvtool_misc.f90 +++ b/src/qcg/solvtool_misc.f90 @@ -299,7 +299,6 @@ subroutine ensemble_lmo(env, fname, self, NTMP, TMPdir, conv) write (jobcall, '(a,1x,a,1x,a,'' --sp --lmo --chrg '',f4.1,1x,a,'' >xtb_lmo.out'')') & & trim(env%ProgName), trim(fname), trim(env%lmover), self%chrg, trim(pipe) k = 0 !counting the finished jobs - !___________________________________________________________________________________ !$omp parallel & @@ -717,7 +716,10 @@ subroutine ens_sp(env, fname, NTMP, TMPdir) do i = 1, NTMP vz = i !$omp task firstprivate( vz ) private( tmppath ) - write (tmppath, '(a,i0)') trim(TMPdir), i + call initsignal() + !$omp critical + write (tmppath, '(a,i0)') trim(TMPdir), vz + !$omp end critical call command('cd '//trim(tmppath)//' && '//trim(jobcall)) !$omp critical k = k + 1