From d7328500376b965f1ef381702cfb09bed94dbc92 Mon Sep 17 00:00:00 2001 From: Philipp Pracht Date: Wed, 17 Jul 2024 00:16:45 +0200 Subject: [PATCH 1/5] Reimplementation of a simple protonation protocol with new calculator infrastructure Signed-off-by: Philipp Pracht --- src/algos/CMakeLists.txt | 1 + src/algos/meson.build | 1 + src/algos/parallel.f90 | 2 + src/algos/protonate.f90 | 210 ++++ src/calculator/api_helpers.F90 | 4 + src/calculator/calc_type.f90 | 53 + src/calculator/gfn0_api.F90 | 27 + src/classes.f90 | 3 +- src/confparse.f90 | 87 +- src/crest_main.f90 | 21 +- src/legacy_algos/deprotonate.f90 | 8 +- src/legacy_algos/pka.f90 | 48 +- src/legacy_algos/protonate.f90 | 16 +- src/legacy_algos/relaxensemble.f90 | 2 +- src/legacy_algos/tautomerize.f90 | 28 +- src/legacy_wrappers.f90 | 13 + src/parsing/parse_calcdata.f90 | 2 + src/propcalc.f90 | 2 +- src/qmhelpers/CMakeLists.txt | 3 + src/qmhelpers/intpack.f90 | 1522 ++++++++++++++++++++++++++++ src/qmhelpers/local.f90 | 1248 +++++++++++++++++++++++ src/qmhelpers/lopt.f90 | 233 +++++ src/qmhelpers/meson.build | 3 + src/sortens.f90 | 2 +- src/strucreader.f90 | 43 + 25 files changed, 3476 insertions(+), 106 deletions(-) create mode 100644 src/algos/protonate.f90 create mode 100644 src/qmhelpers/intpack.f90 create mode 100644 src/qmhelpers/local.f90 create mode 100644 src/qmhelpers/lopt.f90 diff --git a/src/algos/CMakeLists.txt b/src/algos/CMakeLists.txt index 3c0310e7..3df8c266 100644 --- a/src/algos/CMakeLists.txt +++ b/src/algos/CMakeLists.txt @@ -28,6 +28,7 @@ list(APPEND srcs "${dir}/search_1.f90" "${dir}/search_mecp.f90" "${dir}/setuptest.f90" + "${dir}/protonate.f90" "${dir}/hessian_tools.f90" "${dir}/ConfSolv.F90" "${dir}/search_conformers.f90" diff --git a/src/algos/meson.build b/src/algos/meson.build index 48834277..43977333 100644 --- a/src/algos/meson.build +++ b/src/algos/meson.build @@ -26,6 +26,7 @@ srcs += files( 'search_1.f90', 'search_mecp.f90', 'setuptest.f90', + 'protonate.f90', 'hessian_tools.f90', 'ConfSolv.F90', 'search_conformers.f90', diff --git a/src/algos/parallel.f90 b/src/algos/parallel.f90 index 9d0e79ae..45b1e10f 100644 --- a/src/algos/parallel.f90 +++ b/src/algos/parallel.f90 @@ -178,6 +178,8 @@ subroutine crest_oloop(env,nat,nall,at,xyz,eread,dump) !* subroutine crest_oloop !* This subroutine performs concurrent geometry optimizations !* for the given ensemble. Inputs xyz and eread are overwritten +!* +!* IMPORTANT: xyz should be in Bohr(!) for this routine !*************************************************************** use crest_parameters,only:wp,stdout,sep use crest_calculator diff --git a/src/algos/protonate.f90 b/src/algos/protonate.f90 new file mode 100644 index 00000000..26cc0e4b --- /dev/null +++ b/src/algos/protonate.f90 @@ -0,0 +1,210 @@ +!================================================================================! +! This file is part of crest. +! +! Copyright (C) 2022 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 +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! crest is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with crest. If not, see . +!================================================================================! + +!========================================================================================! +!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>! +!> New implementation of protonation routines +!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>! +!========================================================================================! +subroutine crest_new_protonate(env,tim) +!******************************************************************** +!* Standalone runtype for protonating a molecule +!* +!* Input/Output: +!* env - crest's systemdata object +!* tim - timer object +!******************************************************************** + use crest_parameters + use crest_data + use crest_calculator + use strucrd + use optimize_module + implicit none + type(systemdata),intent(inout) :: env + type(timer),intent(inout) :: tim + type(coord) :: mol,molnew + integer :: i,j,k,l,io,ich,T,Tn,np + logical :: pr,wr + character(len=80) :: atmp +!========================================================================================! + type(calcdata) :: calc + real(wp) :: accuracy,etemp + + real(wp) :: energy,dip + real(wp),allocatable :: grad(:,:) + type(calcdata),allocatable :: tmpcalc + real(wp),allocatable :: protxyz(:,:) + integer :: natp + integer,allocatable :: atp(:) + real(wp),allocatable :: xyzp(:,:,:) + real(wp),allocatable :: ep(:) + character(len=*),parameter :: partial = '∂E/∂' +!========================================================================================! + write (stdout,*) + !call system('figlet singlepoint') + write (stdout,*) " _ _ " + write (stdout,*) " _ __ _ __ ___ | |_ ___ _ __ __ _| |_ ___ " + write (stdout,*) "| '_ \| '__/ _ \| __/ _ \| '_ \ / _` | __/ _ \ " + write (stdout,*) "| |_) | | | (_) | || (_) | | | | (_| | || __/ " + write (stdout,*) "| .__/|_| \___/ \__\___/|_| |_|\__,_|\__\___| " + write (stdout,*) "|_| " + write (stdout,*) "-----------------------------------------------" + write (stdout,*) " automated protonation site screening script " + write (stdout,*) 'Cite as:' + write (stdout,*) 'P.Pracht, C.A.Bauer, S.Grimme' + write (stdout,*) 'JCC, 2017, 38, 2618–2631.' + write (stdout,*) + +!========================================================================================! + call new_ompautoset(env,'max',0,T,Tn) + call ompprint_intern() +!========================================================================================! + call env%ref%to(mol) + write (stdout,*) + write (stdout,*) 'Input structure:' + call mol%append(stdout) + write (stdout,*) +!========================================================================================! +!>--- The first step is to perfom a LMO calculation to identify suitable protonation sites + call tim%start(14,'LMO center calculation') + write (stdout,'(a)') repeat('-',80) + write (stdout,'(a)') + write (stdout,'(a)',advance='no') '> Setting up GFN0-xTB for LMO center calculation ... ' + flush (stdout) + + allocate (grad(3,mol%nat),source=0.0_wp) + allocate (tmpcalc) + call env2calc(env,tmpcalc,mol) + tmpcalc%calcs(1)%id = jobtype%gfn0 + tmpcalc%calcs(1)%rdwbo = .true. + tmpcalc%calcs(1)%getlmocent = .true. + tmpcalc%ncalculations = 1 + write (stdout,'(a)') 'done.' +!>--- and then start it + write (stdout,'(a)',advance='yes') '> Performing singlepoint calculation ... ' + flush (stdout) +!>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<< done.' + write (atmp,'(a)') '> Total wall time for calculations' + call tim%write_timing(stdout,14,trim(atmp),.true.) + write (stdout,'(a)') repeat('-',80) + write (stdout,*) + if (io /= 0) then + write (stdout,*) + write (stdout,*) 'WARNING: Calculation exited with error!' + end if +!>--- check LMO center + np = tmpcalc%calcs(1)%nprot + if (np > 0) then + write (stdout,'(a,i0,a)') '> ',np,' π- or LP-centers identified as protonation candidates' + call move_alloc(tmpcalc%calcs(1)%protxyz,protxyz) + !do i=1,np + ! write(stdout,*) protxyz(1:3,i) + !enddo + else + write (stdout,*) + write (stdout,*) 'WARNING: No suitable protonation sites found!' + write (stdout,*) ' Check if you expect π- or LP-centers for your molecule!' + write (stdout,*) + return + end if + deallocate (tmpcalc) + deallocate (grad) + +!========================================================================================! +!>--- If we reached this point, we have candidate positions for our protonation! + write(stdout,'(a)',advance='yes') '> Generating candidate structures ... ' + flush(stdout) + natp = mol%nat+1 + allocate(atp(natp), source=0) + allocate(xyzp(3,natp,np),ep(np),source=0.0_wp) + call protonation_candidates(env,mol,natp,np,protxyz,atp,xyzp) + write(stdout,'(a)') '> Write protonate_0.xyz with candidates ... ' + call wrensemble('protonate_0.xyz',natp,np,atp,xyzp*autoaa,ep) + write(stdout,'(a)') '> done.' + write(stdout,*) + +!>--- Enforce further constraints, conditions, etc. +! (TODO) + +!>--- Optimize candidates + call smallhead('Protomer Ensemble Optimization') + call tim%start(15,'Ensemble optimization') + call print_opt_data(env%calc,stdout) + call crest_oloop(env,natp,np,atp,xyzp,ep,.true.) + call tim%stop(15) + write(stdout,'(a)') '> Write protonate_1.xyz with optimized structures ... ' + call rename(ensemblefile,'protonate_1.xyz') + + +!========================================================================================! + return +end subroutine crest_new_protonate + +!========================================================================================! +!========================================================================================! + +subroutine protonation_candidates(env,mol,natp,np,protxyz,at,xyz) +!******************************************************** +!* generate protonation/ionization candidate structures +!* The outputs are at and xyz, the latter being in Bohr +!******************************************************** + use crest_data + use crest_parameters + use strucrd,only:coord + implicit none + !> INPUT + type(systemdata),intent(inout) :: env + type(coord),intent(in) :: mol + integer,intent(in) :: natp + integer,intent(in) :: np + real(wp),intent(in) :: protxyz(3,np) + !> OUTPUT + integer,intent(out) :: at(natp) + real(wp),intent(out) :: xyz(3,natp,np) + !> LOCAL + integer :: i,j,k,l + + if (natp .ne. mol%nat+1) then + write (stdout,'(a)') 'WARNING: Inconsistent number of atoms in protonation routine' + error stop + end if + + write(stdout,'(a)') '> Increasing the molecular charge by 1' + call env%calc%increase_charge() + +!>--- Populate + do i = 1,np + do j = 1,mol%nat + xyz(1:3,j,i) = mol%xyz(1:3,j) + at(j) = mol%at(j) + end do + xyz(1:3,natp,i) = protxyz(1:3,i) + at(natp) = 1 + end do +end subroutine protonation_candidates + +!========================================================================================! +!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< other properties logical,allocatable :: getsasa(:) + logical :: getlmocent = .false. + integer :: nprot = 0 + real(wp),allocatable :: protxyz(:,:) !>--- API constructs integer :: tblitelvl = 2 @@ -267,6 +270,8 @@ module calc_type procedure :: active => calc_set_active procedure :: active_restore => calc_set_active_restore procedure :: freezegrad => calculation_freezegrad + procedure :: increase_charge => calculation_increase_charge + procedure :: decrease_charge => calculation_decrease_charge end type calcdata public :: get_dipoles @@ -605,6 +610,54 @@ subroutine calculation_freezegrad(self,grad) end if end subroutine calculation_freezegrad +!=========================================================================================! + + subroutine calculation_increase_charge(self,dchrg) +!****************************************************************** +!* increase the charge of all calculation_settings objects by one +!* or the specified dchrg +!****************************************************************** + implicit none + class(calcdata) :: self + integer,intent(in),optional :: dchrg + integer :: i,j + if(self%ncalculations > 0)then + if(present(dchrg))then + j = dchrg + else + j = 1 + endif + do i=1,self%ncalculations + self%calcs(i)%chrg = self%calcs(i)%chrg + j + enddo + endif + return + end subroutine calculation_increase_charge + +!=========================================================================================! + + subroutine calculation_decrease_charge(self,dchrg) +!****************************************************************** +!* decrease the charge of all calculation_settings objects by one +!* or the specified dchrg +!****************************************************************** + implicit none + class(calcdata) :: self + integer,intent(in),optional :: dchrg + integer :: i,j + if(self%ncalculations > 0)then + if(present(dchrg))then + j = dchrg + else + j = 1 + endif + do i=1,self%ncalculations + self%calcs(i)%chrg = self%calcs(i)%chrg - j + enddo + endif + return + end subroutine calculation_decrease_charge + !=========================================================================================! subroutine calculation_settings_addconfig(self,config) diff --git a/src/calculator/gfn0_api.F90 b/src/calculator/gfn0_api.F90 index fef1241a..f683975e 100644 --- a/src/calculator/gfn0_api.F90 +++ b/src/calculator/gfn0_api.F90 @@ -30,6 +30,7 @@ module gfn0_api use gfn0_module,only:gfn0_gbsa_init,generate_config #endif use wiberg_mayer,only:get_wbo,get_wbo_rhf,density_matrix + use mo_localize,only:local implicit none private @@ -50,6 +51,7 @@ module gfn0_api public :: gfn0_gen_occ public :: gfn0_print public :: gfn0_getdipole,gfn0_getqat + public :: gfn0_getlmocent !========================================================================================! !========================================================================================! @@ -285,6 +287,31 @@ subroutine gfn0_getqat(g0calc,mol,qat) #endif end subroutine gfn0_getqat +!========================================================================================! + + subroutine gfn0_getlmocent(g0calc,mol,nprot,protxyz) +!******************************************** +!* Localize the MOs to obtain their centers +!******************************************** + implicit none + type(gfn0_data),intent(inout) :: g0calc + type(coord),intent(in) :: mol + real(wp),allocatable :: z(:) + integer,intent(out) :: nprot + real(wp),intent(out),allocatable :: protxyz(:,:) + nprot=0 +#ifdef WITH_GFN0 + !allocate(z(mol%nat), source=0.0_wp) + call mol%get_z(z) + call local(mol%nat,mol%at,g0calc%basis%nbf,g0calc%basis%nao,g0calc%wfn%ihomo,mol%xyz,z, & + & g0calc%wfn%focc,g0calc%wfn%s,g0calc%wfn%p,g0calc%wfn%C,g0calc%wfn%emo, & + & g0calc%basis%aoat2,g0calc%basis%aoat,g0calc%basis%nprim, & + & g0calc%basis%alp,g0calc%basis%lao,g0calc%basis%cont, & + & nprot,protxyz) +#endif + end subroutine gfn0_getlmocent + + !========================================================================================! !========================================================================================! end module gfn0_api diff --git a/src/classes.f90 b/src/classes.f90 index 337fcb43..666d96dc 100644 --- a/src/classes.f90 +++ b/src/classes.f90 @@ -71,6 +71,7 @@ module crest_data integer,parameter,public :: crest_msreac = 9 integer,parameter,public :: crest_pka = 14 integer,parameter,public :: crest_solv = 15 + integer,parameter,public :: crest_protonate = 16 !>> runtypes with IDs between use non-legacy routines <--- property data objects - type(protobj) :: ptb + type(protobj) :: protb !>--- saved constraints type(constra) :: cts diff --git a/src/confparse.f90 b/src/confparse.f90 index 92d9a9d0..f60f6bd2 100644 --- a/src/confparse.f90 +++ b/src/confparse.f90 @@ -268,17 +268,17 @@ subroutine parseflags(env,arg,nra) env%fixfile = 'none selected' !>--- options for possible property calculations, mainly protonation/deprotonation/taut. tool - env%ptb%popthr = 0.01_wp !> = 1% population - env%ptb%ewin = 30.0_wp !> 30 kcal for protonation - env%ptb%swat = 0 - env%ptb%swchrg = 0 - env%ptb%iter = 1 !> number of iteration cycles for tautomerization - env%ptb%swelem = .false. !> replace H⁺ in protonation routine by something else? - env%ptb%allowFrag = .false. !> allow dissociated Structures? - env%ptb%threshsort = .false. !> use ewin threshold window - env%ptb%protdeprot = .false. !> (tautomerize) do first protonation and then deprotonation - env%ptb%deprotprot = .false. !> (tautomerize) do first deprotonation and then protonation - env%ptb%strictPDT = .false. !> strict mode (i.e. bond constraints) for (de)protonation,tautomerization + env%protb%popthr = 0.01_wp !> = 1% population + env%protb%ewin = 30.0_wp !> 30 kcal for protonation + env%protb%swat = 0 + env%protb%swchrg = 0 + env%protb%iter = 1 !> number of iteration cycles for tautomerization + env%protb%swelem = .false. !> replace H⁺ in protonation routine by something else? + env%protb%allowFrag = .false. !> allow dissociated Structures? + env%protb%threshsort = .false. !> use ewin threshold window + env%protb%protdeprot = .false. !> (tautomerize) do first protonation and then deprotonation + env%protb%deprotprot = .false. !> (tautomerize) do first deprotonation and then protonation + env%protb%strictPDT = .false. !> strict mode (i.e. bond constraints) for (de)protonation,tautomerization env%pclean = .false. !> cleanup option for property mode !>--- options for principal component analysis (PCA) and clustering @@ -430,12 +430,12 @@ subroutine parseflags(env,arg,nra) env%performCross = .false. !skip the genetic crossing env%trackorigin = .false. env%Maxrestart = 1 - env%ptb%ewin = 15.0d0 + env%protb%ewin = 15.0d0 env%gbsa = .true. env%solv = '--alpb h2o' - env%ptb%h_acidic = 0 - call pka_argparse(arg(i+1),env%ptb%h_acidic) - if (env%ptb%h_acidic == -2) env%ptb%pka_baseinp = trim(arg(i+1)) + env%protb%h_acidic = 0 + call pka_argparse(arg(i+1),env%protb%h_acidic) + if (env%protb%h_acidic == -2) env%protb%pka_baseinp = trim(arg(i+1)) case ('-compare') !> flag for comparing two ensembles, analysis tool env%compareens = .true. @@ -461,6 +461,7 @@ subroutine parseflags(env,arg,nra) case ('-protonate') !> protonation tool env%properties = p_protonate + env%crestver = crest_protonate write (*,'(2x,a,'' : automated protonation script'')') trim(arg(i)) exit @@ -496,7 +497,7 @@ subroutine parseflags(env,arg,nra) call xyz2coord(env%ensemblename,'coord') !> write coord from lowest structure env%inputcoords = env%ensemblename !> just for a printout if (argument == '-forall') then - env%ptb%alldivers = .true. + env%protb%alldivers = .true. end if exit @@ -652,16 +653,16 @@ subroutine parseflags(env,arg,nra) case ('-exlig','-exligand','-exchligand') env%properties = p_ligand - env%ptb%infile = trim(arg(1)) + env%protb%infile = trim(arg(1)) ctmp = trim(arg(i+1)) - env%ptb%newligand = trim(ctmp) + env%protb%newligand = trim(ctmp) read (arg(i+2),*,iostat=io) j if (io == 0) then - env%ptb%centeratom = j + env%protb%centeratom = j end if read (arg(i+3),*,iostat=io) j if (io == 0) then - env%ptb%ligand = j + env%protb%ligand = j end if exit @@ -669,18 +670,18 @@ subroutine parseflags(env,arg,nra) !> crest --ab --chrg env%properties = p_acidbase if (index(arg(i),'prep') .ne. 0) then - call pka_argparse2(env,arg(i+1),arg(i+2),env%ptb%pka_mode) + call pka_argparse2(env,arg(i+1),arg(i+2),env%protb%pka_mode) else ctmp = trim(arg(i+1)) inquire (file=ctmp,exist=ex) if (ex) then - env%ptb%pka_acidensemble = trim(ctmp) + env%protb%pka_acidensemble = trim(ctmp) write (*,'(1x,a,a)') 'File used for the acid: ',trim(ctmp) end if ctmp = trim(arg(i+2)) inquire (file=ctmp,exist=ex) if (ex) then - env%ptb%pka_baseensemble = trim(ctmp) + env%protb%pka_baseensemble = trim(ctmp) write (*,'(1x,a,a)') 'File used for the base: ',trim(ctmp) end if end if @@ -1487,7 +1488,7 @@ subroutine parseflags(env,arg,nra) call readl(arg(i+1),xx,j) env%ewin = abs(xx(1)) if (any((/p_protonate,p_deprotonate,p_tautomerize/) == env%properties)) then - env%ptb%ewin = abs(xx(1)) + env%protb%ewin = abs(xx(1)) end if write (*,'(2x,a,1x,a)') trim(arg(i)),trim(arg(i+1)) case ('-rthr') !> set RMSD thr @@ -1581,19 +1582,19 @@ subroutine parseflags(env,arg,nra) case ('-protonate') !> protonation tool env%properties = p_protonate env%autozsort = .false. - env%ptb%threshsort = .true. + env%protb%threshsort = .true. case ('-swel') !> switch out H+ to something else in protonation script if (env%properties .eq. -3) then - call swparse(arg(i+1),env%ptb) + call swparse(arg(i+1),env%protb) end if case ('-deprotonate') !> deprotonation tool env%properties = p_deprotonate env%autozsort = .false. - env%ptb%threshsort = .true. + env%protb%threshsort = .true. case ('-tautomerize') !> tautomerization tool env%properties = p_tautomerize env%autozsort = .false. - env%ptb%threshsort = .true. + env%protb%threshsort = .true. case ('-tautomerize2','-exttautomerize') if (env%properties == p_propcalc) then env%properties = p_tautomerize2 @@ -1601,7 +1602,7 @@ subroutine parseflags(env,arg,nra) call env%addjob(abs(p_tautomerize2)) end if env%autozsort = .false. - env%ptb%threshsort = .true. + env%protb%threshsort = .true. env%runver = 33 env%relax = .true. env%performCross = .false. !> skip the genetic crossing @@ -1614,34 +1615,34 @@ subroutine parseflags(env,arg,nra) env%trackorigin = .false. env%Maxrestart = 1 case ('-trev','-tdp') - env%ptb%deprotprot = .true. !> switch to deprotonation-first mode in tautomerization + env%protb%deprotprot = .true. !> switch to deprotonation-first mode in tautomerization case ('-iter') !> number of Protonation/Deprotonation cycles in Tautomerization call readl(arg(i+1),xx,j) - env%ptb%iter = nint(xx(1)) + env%protb%iter = nint(xx(1)) case ('-texcl','-blacklist') ctmp = trim(arg(i+1)) case ('-strict') - env%ptb%strictPDT = .true. + env%protb%strictPDT = .true. case ('-verystrict','-vstrict') - env%ptb%strictPDT = .false. - env%ptb%fixPDT = .true. + env%protb%strictPDT = .false. + env%protb%fixPDT = .true. case ('-fstrict') - env%ptb%strictPDT = .true. - env%ptb%fixPDT = .true. + env%protb%strictPDT = .true. + env%protb%fixPDT = .true. case ('-corr','-abcorr') - env%ptb%strictPDT = .true. - env%ptb%fixPDT = .true. - env%ptb%ABcorrection = .true. + env%protb%strictPDT = .true. + env%protb%fixPDT = .true. + env%protb%ABcorrection = .true. case ('-pkaensemble') env%preopt = .false. env%presp = .false. - call pka_argparse2(env,arg(i+1),arg(i+2),env%ptb%pka_mode) + call pka_argparse2(env,arg(i+1),arg(i+2),env%protb%pka_mode) case ('-pkaparam') - env%ptb%rdcfer = .true. + env%protb%rdcfer = .true. if (i+1 .le. nra) then ctmp = trim(arg(i+1)) if (ctmp(1:1) .ne. '-') then - env%ptb%cferfile = ctmp + env%protb%cferfile = ctmp end if end if !========================================================================================! @@ -2470,7 +2471,7 @@ subroutine inputcoords(env,arg) !>-- for protonation/deprotonation applications get ref. number of fragments !>-- also get some other structure based info call simpletopo_file('coord',zmol,.false.,.false.,'') - env%ptb%nfrag = zmol%nfrag + env%protb%nfrag = zmol%nfrag call zmol%deallocate() return diff --git a/src/crest_main.f90 b/src/crest_main.f90 index b0961ba8..b5446ce7 100644 --- a/src/crest_main.f90 +++ b/src/crest_main.f90 @@ -140,10 +140,10 @@ program CREST case (p_compare) call compare_ensembles(env) call propquit(tim) - !>--- protonation tool - case (p_protonate) - call protonate(env,tim) - call propquit(tim) +! !>--- protonation tool +! case (p_protonate) +! call protonate(env,tim) +! call propquit(tim) !>--- deprotonation case (p_deprotonate) call deprotonate(env,tim) @@ -188,19 +188,19 @@ program CREST !>--- calculate potential correction for acid/base reaction case (p_acidbase) call tim%start(4,'acid/base') - if (env%ptb%pka_mode == 0) then - call acidbase(env,env%ptb%pka_acidensemble,env%ptb%pka_baseensemble,env%chrg,.true., & + if (env%protb%pka_mode == 0) then + call acidbase(env,env%protb%pka_acidensemble,env%protb%pka_baseensemble,env%chrg,.true., & & .false.,dumfloat,.false.,d3,d4,d5,d6,d7,d8) else - call rewrite_AB_ensemble(env,env%ptb%pka_acidensemble,env%ptb%pka_baseensemble) + call rewrite_AB_ensemble(env,env%protb%pka_acidensemble,env%protb%pka_baseensemble) end if call tim%stop(4) call propquit(tim) !>--- calculate potential correction for acid/base reaction case (p_ligand) call tim%start(4,'') - call ligandtool(env%ptb%infile,env%ptb%newligand, & - & env%ptb%centeratom,env%ptb%ligand) + call ligandtool(env%protb%infile,env%protb%newligand, & + & env%protb%centeratom,env%protb%ligand) call tim%stop(4) call propquit(tim) !>--- wrapper for the thermo routine @@ -290,6 +290,9 @@ program CREST case (crest_ensemblesp) !> singlepoints along ensemble call crest_ensemble_singlepoints(env,tim) + case(crest_protonate) + call protonate(env,tim) + case (crest_test) call crest_playground(env,tim) diff --git a/src/legacy_algos/deprotonate.f90 b/src/legacy_algos/deprotonate.f90 index 19ea909d..47d699bd 100644 --- a/src/legacy_algos/deprotonate.f90 +++ b/src/legacy_algos/deprotonate.f90 @@ -64,10 +64,10 @@ subroutine deprotonate(env,tim) call deprotclean call deprothead - if(.not.allocated(env%ptb%atmap))allocate(env%ptb%atmap(env%nat)) - if(.not.env%ptb%strictPDT .and. .not.env%ptb%fixPDT)then + if(.not.allocated(env%protb%atmap))allocate(env%protb%atmap(env%nat)) + if(.not.env%protb%strictPDT .and. .not.env%protb%fixPDT)then !--- sort the input file (H atoms to the bottom) - call htothebottom('coord',env%chrg,env%nat,env%ptb%atmap) + call htothebottom('coord',env%chrg,env%nat,env%protb%atmap) else !--- or sort AND apply heavy atom bond constraints call PDT_constraints(env) @@ -76,7 +76,7 @@ subroutine deprotonate(env,tim) !--- get some settings call getcwd(thispath) deprotname='deprotonate_0.xyz' - deprot=env%ptb + deprot=env%protb refchrg = env%chrg deprot%newchrg = env%chrg - 1 !increase chrg by one env%chrg = env%chrg - 1 !in the new version all calculations access env%chrg!!! diff --git a/src/legacy_algos/pka.f90 b/src/legacy_algos/pka.f90 index 5652e0e3..76d746d0 100644 --- a/src/legacy_algos/pka.f90 +++ b/src/legacy_algos/pka.f90 @@ -78,7 +78,7 @@ end function pKaPolyFER !---- FER parameter selection parinfo = '' - if (env%ptb%rdCFER) then + if (env%protb%rdCFER) then pkaparam = 'external' else !pkaparam ='none' @@ -109,13 +109,13 @@ end function pKaPolyFER c(3) = -0.11103 ! and checked outlier HClO4 c(4) = 0.0001835 case ('external') - call pka_rdparam(env%ptb%cferfile,nc,c,parinfo) + call pka_rdparam(env%protb%cferfile,nc,c,parinfo) case default !do nothing with the values c = 0.0d0 c(2) = 1.0d0 end select - select case (env%ptb%pka_mode) + select case (env%protb%pka_mode) case (0) !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! env%crestver = 2 env%gbsa = .true. @@ -125,29 +125,29 @@ end function pKaPolyFER if ((env%gfnver .ne. '--gfn2')) then error stop 'Error: pKa tool only available for GFN2/ALPB(H2O) calculations' end if - if (env%ptb%h_acidic == 0) then + if (env%protb%h_acidic == 0) then error stop 'Error: no hydrogen atom selected for acid dissociation' end if !---- very important preperation: resort input so that heavy atoms come first, and then all hydrogens - allocate (env%ptb%atmap(env%nat)) - call htothebottom('coord',env%chrg,env%nat,env%ptb%atmap) + allocate (env%protb%atmap(env%nat)) + call htothebottom('coord',env%chrg,env%nat,env%protb%atmap) !--- check for the selected acidic H atom h = 0 do i = 1,env%nat - if (env%ptb%atmap(i) == env%ptb%h_acidic) h = i + if (env%protb%atmap(i) == env%protb%h_acidic) h = i end do - if ((h .ne. env%ptb%h_acidic).and.(env%ptb%h_acidic > 0)) then + if ((h .ne. env%protb%h_acidic).and.(env%protb%h_acidic > 0)) then write (stdout,'(1x,a,i0,a,i0)') 'The position of the selected acidic hydrogen was updated ', & - & env%ptb%h_acidic,' ---> ',h - env%ptb%h_acidic = h + & env%protb%h_acidic,' ---> ',h + env%protb%h_acidic = h end if !---- first do conformational search in water call largehead2('Calculation of the acid') !--- because the acid must be the input file, it is saved as "coord" call ACID%open('coord') - h = env%ptb%h_acidic + h = env%protb%h_acidic if ((h > 0)) then if (ACID%at(h) .ne. 1) error stop 'selected atom in pKa tool is not a hydrogen atom' end if @@ -173,26 +173,26 @@ end function pKaPolyFER !---- then do the (relaxed) deprotonation call largehead2('Calculation of the base') env%ensemblename = 'none selected' !RESET, IMPORTANT! - env%ptb%threshsort = .true. + env%protb%threshsort = .true. !=================================================================================! !--- write the base input file - if (env%ptb%h_acidic > 0) then !-- for a manually selected H atom + if (env%protb%h_acidic > 0) then !-- for a manually selected H atom BASE%nat = refnat-1 allocate (BASE%at(BASE%nat),source=0) allocate (BASE%xyz(3,BASE%nat),source=0.0_wp) j = 0 do i = 1,ACID%nat - if (i == env%ptb%h_acidic) cycle + if (i == env%protb%h_acidic) cycle j = j+1 BASE%at(j) = ACID%at(i) BASE%xyz(1:3,j) = ACID%xyz(1:3,i) end do - else if (env%ptb%h_acidic == -1) then !-- "automatic" mode + else if (env%protb%h_acidic == -1) then !-- "automatic" mode call deprotonate(env,tim) call rmrfw('deprotonate_') call BASE%open('deprotonated.xyz') - else if (env%ptb%h_acidic == -2) then !-- read-in base file mode - call BASE%open(env%ptb%pka_baseinp) + else if (env%protb%h_acidic == -2) then !-- read-in base file mode + call BASE%open(env%protb%pka_baseinp) end if call BASE%write('base.xyz') env%chrg = refchrg-1 @@ -238,11 +238,11 @@ end function pKaPolyFER write (stdout,*) 'format and contain the final free energies in solution (in Eh) as' write (stdout,*) 'comment line for each structure.' write (stdout,*) - call ACIDENSEMBLE%open(env%ptb%pka_acidensemble) - write (stdout,'(1x,a,a,a,i0,a)') 'Ensemble read for acid: ',env%ptb%pka_acidensemble, & + call ACIDENSEMBLE%open(env%protb%pka_acidensemble) + write (stdout,'(1x,a,a,a,i0,a)') 'Ensemble read for acid: ',env%protb%pka_acidensemble, & & ' with ',ACIDENSEMBLE%nall,' structures' - call BASEENSEMBLE%open(env%ptb%pka_baseensemble) - write (stdout,'(1x,a,a,a,i0,a)') 'Ensemble read for base: ',env%ptb%pka_baseensemble, & + call BASEENSEMBLE%open(env%protb%pka_baseensemble) + write (stdout,'(1x,a,a,a,i0,a)') 'Ensemble read for base: ',env%protb%pka_baseensemble, & & ' with ',BASEENSEMBLE%nall,' structures' write (stdout,'(1x,a)') 'Calculating population average from read-in (free) energies ...' nalla = ACIDENSEMBLE%nall @@ -297,7 +297,7 @@ end function pKaPolyFER dum = GB-GA write (stdout,'(1x,a,f16.8,a,f8.2,a)') 'ΔG =',dum,' Eh,',dum*kcal,' kcal/mol' dG = dum+dE - if (env%ptb%pka_mode .ne. 1) then + if (env%protb%pka_mode .ne. 1) then write (stdout,'(1x,a,f16.8,a,f8.2,a)') 'ΔG+Ecorr =',dG,' Eh,',dG*kcal,' kcal/mol' end if write (stdout,*) @@ -418,8 +418,8 @@ subroutine pka_argparse2(env,str1,str2,h) inquire (file=trim(str2),exist=ex2) if (ex.and.ex2) then h = 1 - env%ptb%pka_acidensemble = trim(str1) - env%ptb%pka_baseensemble = trim(str2) + env%protb%pka_acidensemble = trim(str1) + env%protb%pka_baseensemble = trim(str2) end if return end subroutine pka_argparse2 diff --git a/src/legacy_algos/protonate.f90 b/src/legacy_algos/protonate.f90 index 490f2c6e..864db6f8 100644 --- a/src/legacy_algos/protonate.f90 +++ b/src/legacy_algos/protonate.f90 @@ -41,7 +41,7 @@ end subroutine prothead !-------------------------------------------------------------------------------------------- ! Protonation workflow with GFNn-xTB !-------------------------------------------------------------------------------------------- -subroutine protonate(env,tim) +subroutine protonate_legacy(env,tim) use crest_parameters use crest_data use iomod @@ -80,10 +80,10 @@ end subroutine xtblmo error stop 'protonate unavailable with GFN-FF -> need LMOs' endif - if (.not.allocated(env%ptb%atmap)) allocate (env%ptb%atmap(env%nat)) - if (.not.env%ptb%strictPDT.and..not.env%ptb%fixPDT) then + if (.not.allocated(env%protb%atmap)) allocate (env%protb%atmap(env%nat)) + if (.not.env%protb%strictPDT.and..not.env%protb%fixPDT) then !--- sort the input file (H atoms to the bottom) - call htothebottom('coord',env%chrg,env%nat,env%ptb%atmap) + call htothebottom('coord',env%chrg,env%nat,env%protb%atmap) else !--- or sort AND apply heavy atom bond constraints call PDT_constraints(env) @@ -93,7 +93,7 @@ end subroutine xtblmo call getcwd(thispath) dirn = 'PROT' protname = 'protonate_0.xyz' - prot = env%ptb + prot = env%protb refchrg = env%chrg prot%newchrg = env%chrg+1 !increase chrg by one natp = env%nat+1 !additional proton, Nat is increased by one @@ -205,7 +205,7 @@ end subroutine xtblmo close (ich) end if env%nat = natp-1 !reset nat -end subroutine protonate +end subroutine protonate_legacy !-------------------------------------------------------------------------------------------- ! A quick single point xtb calculation and calculate LMOs @@ -269,7 +269,7 @@ subroutine swelem(iname,env) real(wp),allocatable :: eread(:) integer,allocatable :: at(:) - prot = env%ptb + prot = env%protb nchrg = env%chrg+prot%swchrg prot%newchrg = nchrg @@ -291,7 +291,7 @@ subroutine swelem(iname,env) close (ich) deallocate (at,eread,xyz) - env%ptb = prot + env%protb = prot return end subroutine swelem diff --git a/src/legacy_algos/relaxensemble.f90 b/src/legacy_algos/relaxensemble.f90 index 8cca12ee..25bb7c17 100644 --- a/src/legacy_algos/relaxensemble.f90 +++ b/src/legacy_algos/relaxensemble.f90 @@ -52,7 +52,7 @@ subroutine relaxensemble(fname,env,tim) character(len=128) :: atmp !--- deactivate possible constraints of the strict modi - if(env%ptb%strictPDT .or. env%ptb%fixPDT)then + if(env%protb%strictPDT .or. env%protb%fixPDT)then env%cts%used=.false. endif diff --git a/src/legacy_algos/tautomerize.f90 b/src/legacy_algos/tautomerize.f90 index cd91ee81..a764794e 100644 --- a/src/legacy_algos/tautomerize.f90 +++ b/src/legacy_algos/tautomerize.f90 @@ -67,10 +67,10 @@ subroutine tautomerize(env,tim) error stop 'tautomerize unavailable with GFN-FF -> need LMOs' endif - if (.not.allocated(env%ptb%atmap)) allocate (env%ptb%atmap(env%nat)) - if (.not.env%ptb%strictPDT.and..not.env%ptb%fixPDT) then + if (.not.allocated(env%protb%atmap)) allocate (env%protb%atmap(env%nat)) + if (.not.env%protb%strictPDT.and..not.env%protb%fixPDT) then !--- sort the input file (H atoms to the bottom) - call htothebottom('coord',env%chrg,env%nat,env%ptb%atmap) + call htothebottom('coord',env%chrg,env%nat,env%protb%atmap) else !--- or sort AND apply heavy atom bond constraints call PDT_constraints(env) @@ -80,7 +80,7 @@ subroutine tautomerize(env,tim) call getcwd(thispath) dirn = 'PROT' tautname = 'tautomerize_0.xyz' - taut = env%ptb + taut = env%protb natp = env%nat !backup Nat refchrg = env%chrg @@ -608,9 +608,9 @@ subroutine tautomerize_blacklist(env,fname,nat,atlist) return end if - if (.not.allocated(env%ptb%blacklist)) then - allocate (env%ptb%blacklist(nat)) - env%ptb%blacklist = .false. !none of the atoms is initially blacklisted + if (.not.allocated(env%protb%blacklist)) then + allocate (env%protb%blacklist(nat)) + env%protb%blacklist = .false. !none of the atoms is initially blacklisted end if allocate (xyz(3,nat),dist(nat,nat)) @@ -633,9 +633,9 @@ subroutine tautomerize_blacklist(env,fname,nat,atlist) if (unconstrained(i) == 1) then if (at(i) == 1) then k = minloc(dist(i,:),1) - env%ptb%blacklist(k) = .true. + env%protb%blacklist(k) = .true. else - env%ptb%blacklist(i) = .true. + env%protb%blacklist(i) = .true. end if end if end do @@ -708,7 +708,7 @@ subroutine tautomerize_ext(ensemb,env,tim) call getcwd(thispath) dirn = 'PROT' tautname = 'tautomerize_0.xyz' - taut = env%ptb + taut = env%protb natp = env%nat !backup Nat refchrg = env%chrg @@ -899,10 +899,10 @@ subroutine PDT_constraints(env) integer :: i,h,nh !-- sort all H atoms to the end of the file - call htothebottom('coord',env%chrg,env%nat,env%ptb%atmap) + call htothebottom('coord',env%chrg,env%nat,env%protb%atmap) !-- decide between atom fixing and bond constraints - if (env%ptb%strictPDT.or.env%ptb%fixPDT) then + if (env%protb%strictPDT.or.env%protb%fixPDT) then !-- if -for some reason- there is a constraint already specified, remove it if (env%cts%used) then @@ -911,7 +911,7 @@ subroutine PDT_constraints(env) env%cts%used = .true. !-- if heavy-atom bond constraints shall be specified - if (.not.env%ptb%fixPDT) then + if (.not.env%protb%fixPDT) then write (*,*) 'Strict mode active. Heavy atom bond constraints will be applied.' write (*,'(1x,a,f8.4,a)') 'Selected force constant:',env%forceconst,' Eh' @@ -940,7 +940,7 @@ subroutine PDT_constraints(env) h = zmol%hydrogen() !-- count hydrogen nh = zmol%nat-h call fix_first_X_atoms(nh,env%forceconst,'fixpositions') - if (env%ptb%strictPDT) then + if (env%protb%strictPDT) then call getbmat(zmol,3,env%forceconst) call appendto('bondlengths','fixpositions') end if diff --git a/src/legacy_wrappers.f90 b/src/legacy_wrappers.f90 index d934d6ed..83d88899 100644 --- a/src/legacy_wrappers.f90 +++ b/src/legacy_wrappers.f90 @@ -349,4 +349,17 @@ subroutine trialOPT(env) end subroutine trialOPT !================================================================================! +subroutine protonate(env,tim) + use iso_fortran_env,only:wp => real64 + use crest_data + implicit none + type(systemdata) :: env + type(timer) :: tim + if (env%legacy) then + call protonate_legacy(env,tim) + else + call crest_new_protonate(env,tim) + end if +end subroutine protonate + diff --git a/src/parsing/parse_calcdata.f90 b/src/parsing/parse_calcdata.f90 index 7a57172b..977ae89b 100644 --- a/src/parsing/parse_calcdata.f90 +++ b/src/parsing/parse_calcdata.f90 @@ -447,6 +447,8 @@ subroutine parse_setting_bool(job,key,val) job%rdgrad = val case ('refresh') job%apiclean = val + case ('lmo','lmocent') + job%getlmocent = val case ('print') job%pr = val if (val) job%prch = 999 !> the actual ID will be generated automatically diff --git a/src/propcalc.f90 b/src/propcalc.f90 index 92cfc87f..0bd621e7 100644 --- a/src/propcalc.f90 +++ b/src/propcalc.f90 @@ -37,7 +37,7 @@ subroutine protreffrag(env) !------ get number of fragments for original structure allocate (xyz(3,nat),at(nat),molvec(nat)) call rdcoord('coord',nat,at,xyz) - call mrec(env%ptb%nfrag,xyz,nat,at,molvec) !requires xyz in bohr + call mrec(env%protb%nfrag,xyz,nat,at,molvec) !requires xyz in bohr deallocate (molvec,at,xyz) end associate end subroutine protreffrag diff --git a/src/qmhelpers/CMakeLists.txt b/src/qmhelpers/CMakeLists.txt index c9e40274..76954b44 100644 --- a/src/qmhelpers/CMakeLists.txt +++ b/src/qmhelpers/CMakeLists.txt @@ -18,6 +18,9 @@ set(dir "${CMAKE_CURRENT_SOURCE_DIR}") list(APPEND srcs "${dir}/wiberg_mayer.f90" + "${dir}/intpack.f90" + "${dir}/lopt.f90" + "${dir}/local.f90" ) set(srcs ${srcs} PARENT_SCOPE) diff --git a/src/qmhelpers/intpack.f90 b/src/qmhelpers/intpack.f90 new file mode 100644 index 00000000..85231238 --- /dev/null +++ b/src/qmhelpers/intpack.f90 @@ -0,0 +1,1522 @@ +! This file is part of xtb. +! +! Copyright (C) 2017-2020 Stefan Grimme +! +! xtb is free software: you can redistribute it and/or modify it under +! the terms of the GNU Lesser General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! xtb is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with xtb. If not, see . + +!! ------------------------------------------------------------------------ +module xtb_intpack +! use xtb_mctc_accuracy, only : wp + use iso_fortran_env, only: wp=>real64 + +contains + +!! --------------------------------------------------------------[SAW1710]- +! this chunk used to have 138 goto marks and 68 goto statements, +! I replaced all of them with normal do loops and got rid of most +! spaceship operaters, if I could figure out what they were +! supposed to do. +! Now down to 30 goto's and 75 goto marks + +subroutine pola(a,b,ga,gb,gm,gm2,iff1,iff2,iall,aa,bb,lmnexp,lmnfak,est,arg,va) + implicit none + ! aufpunkte,intarray + real(wp) :: a(3),b(3),va,est,arg,lmnfak(84) + real(wp) :: gm,gm2,ga,gb + integer iall,lmnexp(84),iff1,iff2 + ! local + real(wp) :: e(3),d(3),olap3,efact,val + real(wp) :: aa(20),bb(20),dd(84) + integer i + ! if you want f-functions: dimension aa(20),bb(20),dd(84) + + ! --- a,b are centres of gaussians + ! --- ga,gb are their exponents + ! gama=ga+gb + ! gm =1.0d0/gama + ! gm2=0.5*gm + ! --- apply product theorem + ! ----- e is center of product gaussian with exponent gama + do i=1,3 + e(i)=(ga*a(i)+gb*b(i))*gm + enddo + ! --- calculate cartesian prefactor for first gaussian + call rhftce(aa,a,e,iff1) + ! --- calculate cartesian prefactor for second gaussian + call rhftce(bb,b,e,iff2) + ! --- form their product + ! dd = 0 + call prod(aa,bb,dd,iff1,iff2) + val = 0 + do i=1,iall + olap3=lmnfak(i)*arg*gm2**lmnexp(i) + val=dd(i)*olap3+val + enddo + va=exp(-est)*val + +end subroutine pola + + +!======================================================================= +! cartesian gaussian functions (6d,10f...) +! iff : +! s,px, py pz, dx**2 dy**2 dz**2 dxy dxz dyz +! 1 2 3 4 5 6 7 8 9 10 +! fxxx, fyyy, fzzz, fxxy, fxxz, fyyx, fyyz, fxzz, fyzz, fxyz +! 11 12 13 14 15 16 17 18 19 20 +! a, are aufpunkte +! nt : # of returns +! va : integral +!======================================================================= + +!! --------------------------------------------------------------[SAW1710]- +! changed do loops, replaced spaceships +subroutine prola(aname,a,b,etaij4,etakl4,iff1,iff2,va,nt) + implicit real(wp)(a-h,o-z) + external aname + ! aufpunkte,ref point,intarray + real(wp) a(3),b(3),va(nt) + ! local + dimension d(3),dd(84),v(3),val(3),ra(3),rb(3) + dimension e(3),aa(20),bb(20) + integer,parameter :: lin(84) = & + & (/0,1,0,0,2,0,0,1,1,0,3,0,0,2,2,1,0,1,0,1,4,0,0,3,3,1,0,1,0,2,2,0, & + & 2,1,1,5,0,0,3,3,2,2,0,0,4,4,1,0,0,1,1,3,1,2,2,1,6,0,0,3,3,0,5,5, & + & 1,0,0,1,4,4,2,0,2,0,3,3,1,2,2,1,4,1,1,2/) + integer,parameter :: min(84) = & + & (/0,0,1,0,0,2,0,1,0,1,0,3,0,1,0,2,2,0,1,1,0,4,0,1,0,3,3,0,1,2,0,2, & + & 1,2,1,0,5,0,2,0,3,0,3,2,1,0,4,4,1,0,1,1,3,2,1,2,0,6,0,3,0,3,1,0, & + & 0,1,5,5,2,0,0,2,4,4,2,1,3,1,3,2,1,4,1,2/) + integer,parameter :: nin(84) = & + & (/0,0,0,1,0,0,2,0,1,1,0,0,3,0,1,0,1,2,2,1,0,0,4,0,1,0,1,3,3,0,2,2, & + & 1,1,2,0,0,5,0,2,0,3,2,3,0,1,0,1,4,4,3,1,1,1,2,2,0,0,6,0,3,3,0,1, & + & 5,5,1,0,0,2,4,4,0,2,1,2,2,3,1,3,1,1,4,2/) + ! --- a,b are centres of gaussians + ra = a + rb = b + ! --- ga,gb are their exponents + ga=etaij4 + gb=etakl4 + aa = 0 + bb = 0 + aa(iff1)=1.0d0 + bb(iff2)=1.0d0 + ia=iff1 + ib=iff2 + ! --- ia,ib are the canonical indices for monom + ! --- apply product theorem + cij=0.0d0 + ckl=0.0d0 + call divpt(a,etaij4,b,etakl4,cij,ckl,e,gama,efact) + ! if(mprp.eq.16) goto 200 + + ! --- calculate cartesian prefactor for first gaussian + call rhftce(aa,a,e,iff1) + ! --- calculate cartesian prefactor for second gaussian + call rhftce(bb,b,e,iff2) + ! --- form their product + call prod(aa,bb,dd,iff1,iff2) + val = 0 + v = 0 + ! ----- e is center of product gaussian with exponent gama + ! ----- c is reference point + ! d = e - c + d = e + ! aname represents an external function + if(iff1.gt.10.or.iff2.gt.10) goto 110 + if(iff1.gt.4.or.iff2.gt.4) goto 120 + if(iff1.gt.1.or.iff2.gt.1) goto 130 + ! s-s - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + call aname(ia,ib,ga,gb,ra,rb,lin(1),min(1),nin(1),gama,v,d) + do j=1,nt + val(j)=dd(1)*v(j)+val(j) + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + 130 continue + if(iff1.ne.1.and.iff2.ne.1) goto 131 + ! s-p - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + do i=1,4 + if(dabs(dd(i)).le.1.d-8) cycle + call aname(ia,ib,ga,gb,ra,rb,lin(i),min(i),nin(i),gama,v,d) + do j=1,nt + val(j)=dd(i)*v(j)+val(j) + enddo + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + + 131 continue + ! p-p - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + do i=1,10 + if(dabs(dd(i)).le.1.d-8) cycle + call aname(ia,ib,ga,gb,ra,rb,lin(i),min(i),nin(i),gama,v,d) + do j=1,nt + val(j)=dd(i)*v(j)+val(j) + enddo + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + + 120 continue + if(iff1.ne.1.and.iff2.ne.1) goto 121 + ! s-d - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + do i=1,10 + if(dabs(dd(i)).le.1.d-8) cycle + call aname(ia,ib,ga,gb,ra,rb,lin(i),min(i),nin(i),gama,v,d) + do j=1,nt + val(j)=dd(i)*v(j)+val(j) + enddo + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + + 121 continue + if(iff1.gt.4.and.iff2.gt.4) goto 122 + ! p-d - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + do i=1,20 + if(dabs(dd(i)).le.1.d-8) cycle + call aname(ia,ib,ga,gb,ra,rb,lin(i),min(i),nin(i),gama,v,d) + do j=1,nt + val(j)=dd(i)*v(j)+val(j) + enddo + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + + 122 continue + ! d-d - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + do i=1,35 + if(dabs(dd(i)).le.1.d-8) cycle + call aname(ia,ib,ga,gb,ra,rb,lin(i),min(i),nin(i),gama,v,d) + do j=1,nt + val(j)=dd(i)*v(j)+val(j) + enddo + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + + 110 continue + if(iff1.ne.1.and.iff2.ne.1) goto 111 + ! s-f - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + do i=1,20 + if(dabs(dd(i)).le.1.d-8) cycle + call aname(ia,ib,ga,gb,ra,rb,lin(i),min(i),nin(i),gama,v,d) + do j=1,nt + val(j)=dd(i)*v(j)+val(j) + enddo + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + + 111 continue + if(iff1.gt.4.and.iff2.gt.4) goto 112 + ! p-f - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + do i=1,35 + if(dabs(dd(i)).le.1.d-8) cycle + call aname(ia,ib,ga,gb,ra,rb,lin(i),min(i),nin(i),gama,v,d) + do j=1,nt + val(j)=dd(i)*v(j)+val(j) + enddo + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + + 112 continue + if(iff1.gt.10.and.iff2.gt.10) goto 113 + ! d-f - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + do i=1,56 + if(dabs(dd(i)).le.1.d-8) cycle + call aname(ia,ib,ga,gb,ra,rb,lin(i),min(i),nin(i),gama,v,d) + do j=1,nt + val(j)=dd(i)*v(j)+val(j) + enddo + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + 113 continue + ! f-f - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + do i=1,84 + if(dabs(dd(i)).le.1.d-8) cycle + call aname(ia,ib,ga,gb,ra,rb,lin(i),min(i),nin(i),gama,v,d) + do j=1,nt + val(j)=dd(i)*v(j)+val(j) + enddo + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + +end subroutine prola + + +!======================================================================= +! cartesian gaussian functions (6d,10f...) +! iff : +! s,px, py pz, dx**2 dy**2 dz**2 dxy dxz dyz +! 1 2 3 4 5 6 7 8 9 10 +! fxxx, fyyy, fzzz, fxxy, fxxz, fyyx, fyyz, fxzz, fyzz, fxyz +! 11 12 13 14 15 16 17 18 19 20 +! a, are aufpunkte +! nt : # of returns +! va : integral +!======================================================================= + +!! --------------------------------------------------------------[SAW1710]- +! changed do loops, replaced spaceships +subroutine propa(aname,a,b,c,etaij4,etakl4,iff1,iff2,va,nt) + implicit real(wp)(a-h,o-z) + external aname + ! aufpunkte,ref point,intarray + real(wp) a(3),b(3),c(3),va(nt) + ! local + common /abfunc/ ra(3),rb(3),ga,gb,ia,ib + dimension d(3),dd(84),v(nt),val(nt) !,v(3),val(3) + dimension e(3),aa(20),bb(20) + integer,parameter :: lin(84) = & + & (/0,1,0,0,2,0,0,1,1,0,3,0,0,2,2,1,0,1,0,1,4,0,0,3,3,1,0,1,0,2,2,0, & + & 2,1,1,5,0,0,3,3,2,2,0,0,4,4,1,0,0,1,1,3,1,2,2,1,6,0,0,3,3,0,5,5, & + & 1,0,0,1,4,4,2,0,2,0,3,3,1,2,2,1,4,1,1,2/) + integer,parameter :: min(84) = & + & (/0,0,1,0,0,2,0,1,0,1,0,3,0,1,0,2,2,0,1,1,0,4,0,1,0,3,3,0,1,2,0,2, & + & 1,2,1,0,5,0,2,0,3,0,3,2,1,0,4,4,1,0,1,1,3,2,1,2,0,6,0,3,0,3,1,0, & + & 0,1,5,5,2,0,0,2,4,4,2,1,3,1,3,2,1,4,1,2/) + integer,parameter :: nin(84) = & + & (/0,0,0,1,0,0,2,0,1,1,0,0,3,0,1,0,1,2,2,1,0,0,4,0,1,0,1,3,3,0,2,2, & + & 1,1,2,0,0,5,0,2,0,3,2,3,0,1,0,1,4,4,3,1,1,1,2,2,0,0,6,0,3,3,0,1, & + & 5,5,1,0,0,2,4,4,0,2,1,2,2,3,1,3,1,1,4,2/) + + ! --- a,b are centres of gaussians + ra = a + rb = b + ! --- ga,gb are their exponents + ga=etaij4 + gb=etakl4 + aa = 0 + bb = 0 + aa(iff1)=1.0d0 + bb(iff2)=1.0d0 + ia=iff1 + ib=iff2 + ! --- ia,ib are the canonical indices for monom + ! --- apply product theorem + cij=0.0d0 + ckl=0.0d0 + call divpt(a,etaij4,b,etakl4,cij,ckl,e,gama,efact) + ! if(mprp.eq.16) goto 200 + + ! --- calculate cartesian prefactor for first gaussian + call rhftce(aa,a,e,iff1) + ! --- calculate cartesian prefactor for second gaussian + call rhftce(bb,b,e,iff2) + ! --- form their product + call prod(aa,bb,dd,iff1,iff2) + val = 0 + ! ----- e is center of product gaussian with exponent gama + ! ----- c is reference point + d = e - c + ! d = e + ! aname represents an external function + if(iff1.gt.10.or.iff2.gt.10) goto 110 + if(iff1.gt.4.or.iff2.gt.4) goto 120 + if(iff1.gt.1.or.iff2.gt.1) goto 130 + ! s-s - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + call aname(lin(1),min(1),nin(1),gama,v,d) + do j=1,nt + val(j)=dd(1)*v(j)+val(j) + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + 130 continue + if(iff1.ne.1.and.iff2.ne.1) goto 131 + ! s-p - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + do i=1,4 + if(dabs(dd(i)).le.1.d-8) cycle + call aname(lin(i),min(i),nin(i),gama,v,d) + do j=1,nt + val(j)=dd(i)*v(j)+val(j) + enddo + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + 131 continue + ! p-p - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + do i=1,10 + if(dabs(dd(i)).le.1.d-8) cycle + call aname(lin(i),min(i),nin(i),gama,v,d) + do j=1,nt + val(j)=dd(i)*v(j)+val(j) + enddo + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + 120 continue + if(iff1.ne.1.and.iff2.ne.1) goto 121 + ! s-d - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + do i=1,10 + if(dabs(dd(i)).le.1.d-8) cycle + call aname(lin(i),min(i),nin(i),gama,v,d) + do j=1,nt + val(j)=dd(i)*v(j)+val(j) + enddo + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + 121 continue + if(iff1.gt.4.and.iff2.gt.4) goto 122 + ! p-d - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + do i=1,20 + if(dabs(dd(i)).le.1.d-8) cycle + call aname(lin(i),min(i),nin(i),gama,v,d) + do j=1,nt + val(j)=dd(i)*v(j)+val(j) + enddo + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + 122 continue + ! d-d - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + do i=1,35 + if(dabs(dd(i)).le.1.d-8) cycle + call aname(lin(i),min(i),nin(i),gama,v,d) + do j=1,nt + val(j)=dd(i)*v(j)+val(j) + enddo + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + 110 continue + if(iff1.ne.1.and.iff2.ne.1) goto 111 + ! s-f - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + do i=1,20 + if(dabs(dd(i)).le.1.d-8) cycle + call aname(lin(i),min(i),nin(i),gama,v,d) + do j=1,nt + val(j)=dd(i)*v(j)+val(j) + enddo + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + 111 continue + if(iff1.gt.4.and.iff2.gt.4) goto 112 + ! p-f - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + do i=1,35 + if(dabs(dd(i)).le.1.d-8) cycle + call aname(lin(i),min(i),nin(i),gama,v,d) + do j=1,nt + val(j)=dd(i)*v(j)+val(j) + enddo + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + 112 continue + if(iff1.gt.10.and.iff2.gt.10) goto 113 + ! d-f - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + do i=1,56 + if(dabs(dd(i)).le.1.d-8) cycle + call aname(lin(i),min(i),nin(i),gama,v,d) + do j=1,nt + val(j)=dd(i)*v(j)+val(j) + enddo + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + 113 continue + ! f-f - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + do i=1,84 + if(dabs(dd(i)).le.1.d-8) cycle + call aname(lin(i),min(i),nin(i),gama,v,d) + do j=1,nt + val(j)=dd(i)*v(j)+val(j) + enddo + enddo + do i=1,nt + va(i)=efact*val(i) + enddo + return + + ! ang mom + !cha + ! 200 do i=1,3 + ! va(i)=0.d0 + ! d(i)=e(i)-c(i) + ! enddo + ! if (efact.eq.0.d0) return + ! call aname(lll,mmm,nnn,gama,v,d) + ! efact=aa(ia)*bb(ib)*efact + ! do i=1,nt + ! va(i)=v(i)*efact + ! enddo + ! return + !cha +end subroutine propa + + +!! --------------------------------------------------------------[SAW1710]- +! changed do loops +pure subroutine divpt(a,alpha,b,beta,ca,cb,e,fgama,ffact) + implicit none + ! this computes the center, exponent, and multiplying factor of + ! a single gaussian which can replace the product of two gaussia + ! centers a and b, and exponents alpha and beta. + real(wp), intent(in) :: a(3),alpha,b(3),beta,ca,cb + real(wp), intent(out) :: fgama,ffact,e(3) + integer :: i + real(wp) :: absqd,abprod,tol + real(wp), parameter :: toluol = 20.7232660d0 + ! abexp=alpha+beta + do i=1,3 + e(i)=(alpha*a(i)+beta*b(i))/(alpha+beta) + enddo + absqd=0.0d0 + do i=1,3 + absqd=absqd+(b(i)-a(i))*(b(i)-a(i)) + enddo + abprod=absqd*alpha*beta/(alpha+beta) + fgama=(alpha+beta) + ffact=0.d0 + tol=toluol+ca+cb + if(abprod.gt.tol) return + ffact=dexp(-abprod) + return +end subroutine divpt + +!! --------------------------------------------------------------[SAW1710]- +! made pure, made explicit +pure subroutine rhftce(cfs,a,e,iff) + implicit none + integer,intent(in) :: iff + real(wp), intent(in) :: a(*),e(*) + real(wp), intent(inout) :: cfs(*) + real(wp), parameter :: c2 = 2.0d0 + real(wp), parameter :: c3 = 3.0d0 + real(wp) :: aex,aey,aez + ! ---- e = center of product function, a = center of single gaussian + aex = e(1)-a(1) + aey = e(2)-a(2) + aez = e(3)-a(3) + + select case(iff) + case(1) + continue + case(2) + cfs(1)=aex*cfs(2) + case(3) + cfs(1)=aey*cfs(3) + case(4) + cfs(1)=aez*cfs(4) + case(5) + cfs(1)=aex*aex*cfs(5) + cfs(2)=c2*aex*cfs(5) + case(6) + cfs(1)=aey*aey*cfs(6) + cfs(3)=c2*aey*cfs(6) + case(7) + cfs(1)=aez*aez*cfs(7) + cfs(4)=c2*aez*cfs(7) + case(8) + cfs(1)=aex*aey*cfs(8) + cfs(2)=aey*cfs(8) + cfs(3)=aex*cfs(8) + case(9) + cfs(1)=aex*aez*cfs(9) + cfs(2)=aez*cfs(9) + cfs(4)=aex*cfs(9) + case(10) + cfs(1)=aey*aez*cfs(10) + cfs(3)=aez*cfs(10) + cfs(4)=aey*cfs(10) + case(11) + cfs(1)=aex*aex*aex*cfs(11) + cfs(2)=c3*aex*aex*cfs(11) + cfs(5)=c3*aex*cfs(11) + case(12) + cfs(1)=aey*aey*aey*cfs(12) + cfs(3)=c3*aey*aey*cfs(12) + cfs(6)=c3*aey*cfs(12) + case(13) + cfs(1)=aez*aez*aez*cfs(13) + cfs(4)=c3*aez*aez*cfs(13) + cfs(7)=c3*aez*cfs(13) + case(14) + cfs(1)=aex*aex*aey*cfs(14) + cfs(2)=c2*aex*aey*cfs(14) + cfs(3)=aex*aex*cfs(14) + cfs(5)=aey*cfs(14) + cfs(8)=c2*aex*cfs(14) + case(15) + cfs(1)=aex*aex*aez*cfs(15) + cfs(2)=c2*aex*aez*cfs(15) + cfs(4)=aex*aex*cfs(15) + cfs(5)=aez*cfs(15) + cfs(9)=c2*aex*cfs(15) + case(16) + cfs(1)=aey*aey*aex*cfs(16) + cfs(2)=aey*aey*cfs(16) + cfs(3)=c2*aey*aex*cfs(16) + cfs(6)=aex*cfs(16) + cfs(8)=c2*aey*cfs(16) + case(17) + cfs(1)=aey*aey*aez*cfs(17) + cfs(3)=c2*aey*aez*cfs(17) + cfs(4)=aey*aey*cfs(17) + cfs(6)=aez*cfs(17) + cfs(10)=c2*aey*cfs(17) + case(18) + cfs(1)=aez*aez*aex*cfs(18) + cfs(2)=aez*aez*cfs(18) + cfs(4)=c2*aez*aex*cfs(18) + cfs(7)=aex*cfs(18) + cfs(9)=c2*aez*cfs(18) + case(19) + cfs(1)=aez*aez*aey*cfs(19) + cfs(3)=aez*aez*cfs(19) + cfs(4)=c2*aez*aey*cfs(19) + cfs(7)=aey*cfs(19) + cfs(10)=c2*aez*cfs(19) + case(20) + cfs(1)=aex*aey*aez*cfs(20) + cfs(2)=aez*aey*cfs(20) + cfs(3)=aex*aez*cfs(20) + cfs(4)=aex*aey*cfs(20) + cfs(8)=aez*cfs(20) + cfs(9)=aey*cfs(20) + cfs(10)=aex*cfs(20) + case default + continue + end select + + return +end subroutine rhftce + + +!! --------------------------------------------------------------[SAW1710]- +! made explicit, made pure +pure subroutine prod(c,d,s,iff1,iff2) + implicit none + integer,intent(in) :: iff1,iff2 + real(wp), intent(in) :: c(*),d(*) + real(wp), intent(out) :: s(*) + if(iff1.gt.10.or.iff2.gt.10) goto 30 + if(iff1.gt.4.or.iff2.gt.4) goto 20 + ! s-s - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + s( 1)=c( 1)*d( 1) + ! end of s - s + if(iff1.eq.1.and.iff2.eq.1) return + ! s-p - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + s( 2)=c( 1)*d( 2)+c( 2)*d( 1) + s( 3)=c( 1)*d( 3)+c( 3)*d( 1) + s( 4)=c( 1)*d( 4)+c( 4)*d( 1) + ! end of s - p + if(iff1.eq.1.or.iff2.eq.1) return + ! p-p - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + s( 5)=c( 2)*d( 2) + s( 6)=c( 3)*d( 3) + s( 7)=c( 4)*d( 4) + s( 8)=c( 2)*d( 3)+c( 3)*d( 2) + s( 9)=c( 2)*d( 4)+c( 4)*d( 2) + s(10)=c( 3)*d( 4)+c( 4)*d( 3) + ! end of p - p + return + 20 continue + ! s-d - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + s( 1)=c( 1)*d( 1) + s( 2)=c( 1)*d( 2)+c( 2)*d( 1) + s( 3)=c( 1)*d( 3)+c( 3)*d( 1) + s( 4)=c( 1)*d( 4)+c( 4)*d( 1) + s( 5)=c( 1)*d( 5)+c( 5)*d( 1)+c( 2)*d( 2) + s( 6)=c( 1)*d( 6)+c( 6)*d( 1)+c( 3)*d( 3) + s( 7)=c( 1)*d( 7)+c( 7)*d( 1)+c( 4)*d( 4) + s( 8)=c( 1)*d( 8)+c( 8)*d( 1)+c( 2)*d( 3)+c( 3)*d( 2) + s( 9)=c( 1)*d( 9)+c( 9)*d( 1)+c( 2)*d( 4)+c( 4)*d( 2) + s(10)=c( 1)*d(10)+c(10)*d( 1)+c( 3)*d( 4)+c( 4)*d( 3) + ! end of s - d + if(iff1.eq.1.or.iff2.eq.1) return + ! p-d - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + s(11)=c( 2)*d( 5)+c( 5)*d( 2) + s(12)=c( 3)*d( 6)+c( 6)*d( 3) + s(13)=c( 4)*d( 7)+c( 7)*d( 4) + s(14)=c( 2)*d( 8)+c( 8)*d( 2)+c( 3)*d( 5)+c( 5)*d( 3) + s(15)=c( 2)*d( 9)+c( 9)*d( 2)+c( 4)*d( 5)+c( 5)*d( 4) + s(16)=c( 2)*d( 6)+c( 6)*d( 2)+c( 3)*d( 8)+c( 8)*d( 3) + s(17)=c( 3)*d(10)+c(10)*d( 3)+c( 4)*d( 6)+c( 6)*d( 4) + s(18)=c( 2)*d( 7)+c( 7)*d( 2)+c( 4)*d( 9)+c( 9)*d( 4) + s(19)=c( 3)*d( 7)+c( 7)*d( 3)+c( 4)*d(10)+c(10)*d( 4) + s(20)=c( 2)*d(10)+c(10)*d( 2)+c( 3)*d( 9) & + & +c( 9)*d( 3)+c( 4)*d( 8)+c( 8)*d( 4) + ! end of p - d + if(iff1.lt.5.or.iff2.lt.5) return + ! d-d - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + s(21)=c( 5)*d( 5) + s(22)=c( 6)*d( 6) + s(23)=c( 7)*d( 7) + s(24)=c( 5)*d( 8)+c( 8)*d( 5) + s(25)=c( 5)*d( 9)+c( 9)*d( 5) + s(26)=c( 6)*d( 8)+c( 8)*d( 6) + s(27)=c( 6)*d(10)+c(10)*d( 6) + s(28)=c( 7)*d( 9)+c( 9)*d( 7) + s(29)=c( 7)*d(10)+c(10)*d( 7) + s(30)=c( 5)*d( 6)+c( 6)*d( 5)+c( 8)*d( 8) + s(31)=c( 5)*d( 7)+c( 7)*d( 5)+c( 9)*d( 9) + s(32)=c( 6)*d( 7)+c( 7)*d( 6)+c(10)*d(10) + s(33)=c( 5)*d(10)+c(10)*d( 5)+c( 8)*d( 9)+c( 9)*d( 8) + s(34)=c( 6)*d( 9)+c( 9)*d( 6)+c( 8)*d(10)+c(10)*d( 8) + s(35)=c( 7)*d( 8)+c( 8)*d( 7)+c( 9)*d(10)+c(10)*d( 9) + ! end of d - d + return + 30 continue + ! s-f - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + s( 1)=c( 1)*d( 1) + s( 2)=c( 1)*d( 2)+c( 2)*d( 1) + s( 3)=c( 1)*d( 3)+c( 3)*d( 1) + s( 4)=c( 1)*d( 4)+c( 4)*d( 1) + s( 5)=c( 1)*d( 5)+c( 5)*d( 1)+c( 2)*d( 2) + s( 6)=c( 1)*d( 6)+c( 6)*d( 1)+c( 3)*d( 3) + s( 7)=c( 1)*d( 7)+c( 7)*d( 1)+c( 4)*d( 4) + s( 8)=c( 1)*d( 8)+c( 8)*d( 1)+c( 2)*d( 3)+c( 3)*d( 2) + s( 9)=c( 1)*d( 9)+c( 9)*d( 1)+c( 2)*d( 4)+c( 4)*d( 2) + s(10)=c( 1)*d(10)+c(10)*d( 1)+c( 3)*d( 4)+c( 4)*d( 3) + s(11)=c( 2)*d( 5)+c( 5)*d( 2)+c(1)*d(11)+c(11)*d(1) + s(12)=c( 3)*d( 6)+c( 6)*d( 3)+c(1)*d(12)+c(12)*d(1) + s(13)=c( 4)*d( 7)+c( 7)*d( 4)+c(1)*d(13)+c(13)*d(1) + s(14)=c( 2)*d( 8)+c( 8)*d( 2)+c( 3)*d( 5)+c( 5)*d( 3)+c(1)*d(14)+ & + & c(14)*d(1) + s(15)=c( 2)*d( 9)+c( 9)*d( 2)+c( 4)*d( 5)+c( 5)*d( 4)+c(1)*d(15)+ & + & c(15)*d(1) + s(16)=c( 2)*d( 6)+c( 6)*d( 2)+c( 3)*d( 8)+c( 8)*d( 3)+c(1)*d(16)+ & + & c(16)*d(1) + s(17)=c( 3)*d(10)+c(10)*d( 3)+c( 4)*d( 6)+c( 6)*d( 4)+c(1)*d(17)+ & + & c(17)*d(1) + s(18)=c( 2)*d( 7)+c( 7)*d( 2)+c( 4)*d( 9)+c( 9)*d( 4)+c(1)*d(18)+ & + & c(18)*d(1) + s(19)=c( 3)*d( 7)+c( 7)*d( 3)+c( 4)*d(10)+c(10)*d( 4)+c(1)*d(19)+ & + & c(19)*d(1) + s(20)=c( 2)*d(10)+c(10)*d( 2)+c( 3)*d( 9)+c(9)*d(3)+c(4)*d(8)+ & + & c(8)*d(4)+c(1)*d(20)+c(20)*d(1) + ! end of s - f + if(iff1.eq.1.or.iff2.eq.1) return + ! p-f - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + s(21)=c( 5)*d( 5)+c(2)*d(11)+c(11)*d(2) + s(22)=c( 6)*d( 6)+c(3)*d(12)+c(12)*d(3) + s(23)=c( 7)*d( 7)+c(4)*d(13)+c(13)*d(4) + s(24)=c( 5)*d( 8)+c( 8)*d( 5)+c(3)*d(11)+c(11)*d(3)+c(2)*d(14)+ & + & c(14)*d(2) + s(25)=c( 5)*d( 9)+c( 9)*d( 5)+c(2)*d(15)+c(15)*d(2)+c(4)*d(11)+ & + & c(11)*d(4) + s(26)=c( 6)*d( 8)+c( 8)*d( 6)+c(2)*d(12)+c(12)*d(2)+c(3)*d(16)+ & + & c(16)*d(3) + s(27)=c( 6)*d(10)+c(10)*d( 6)+c(3)*d(17)+c(17)*d(3)+c(4)*d(12)+ & + & c(12)*d(4) + s(28)=c( 7)*d( 9)+c( 9)*d( 7)+c(2)*d(13)+c(13)*d(2)+c(4)*d(18)+ & + & c(18)*d(4) + s(29)=c( 7)*d(10)+c(10)*d( 7)+c(3)*d(13)+c(13)*d(3)+c(4)*d(19)+ & + & c(19)*d(4) + s(30)=c( 5)*d( 6)+c( 6)*d( 5)+c( 8)*d( 8)+c(2)*d(16)+c(16)*d(2)+ & + & c(3)*d(14)+c(14)*d(3) + s(31)=c( 5)*d( 7)+c( 7)*d( 5)+c( 9)*d( 9)+c(2)*d(18)+c(18)*d(2)+ & + & c(4)*d(15)+c(15)*d(4) + s(32)=c( 6)*d( 7)+c( 7)*d( 6)+c(10)*d(10)+c(3)*d(19)+c(19)*d(3)+ & + & c(4)*d(17)+c(17)*d(4) + s(33)=c( 5)*d(10)+c(10)*d( 5)+c( 8)*d( 9)+c( 9)*d( 8)+c(3)*d(15)+ & + & c(15)*d(3)+c(4)*d(14)+c(14)*d(4)+c(2)*d(20)+c(20)*d(2) + s(34)=c( 6)*d( 9)+c( 9)*d( 6)+c( 8)*d(10)+c(10)*d( 8)+c(2)*d(17)+ & + & d(2)*c(17)+c(3)*d(20)+c(20)*d(3)+c(4)*d(16)+c(16)*d(4) + s(35)=c( 7)*d( 8)+c( 8)*d( 7)+c( 9)*d(10)+c(10)*d( 9)+c(2)*d(19)+ & + & c(19)*d(2)+c(3)*d(18)+c(18)*d(3)+c(4)*d(20)+c(20)*d(4) + ! end of p - f + if(iff1.eq.2.or.iff2.eq.2) return + ! d-f - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + s(36)=c(5)*d(11)+c(11)*d(5) + s(37)=c(6)*d(12)+c(12)*d(6) + s(38)=c(7)*d(13)+c(13)*d(7) + s(39)=c(6)*d(11)+c(11)*d(6)+c(5)*d(16)+c(16)*d(5)+c(8)*d(14)+ & + & c(14)*d(8) + s(40)=c(7)*d(11)+c(11)*d(7)+c(5)*d(18)+c(18)*d(5)+c(9)*d(15)+ & + & c(15)*d(9) + s(41)=c(5)*d(12)+c(12)*d(5)+c(6)*d(14)+c(14)*d(6)+c(8)*d(16)+ & + & c(16)*d(8) + s(42)=c(5)*d(13)+c(13)*d(5)+c(7)*d(15)+c(15)*d(7)+c(9)*d(18)+ & + & c(18)*d(9) + s(43)=c(7)*d(12)+c(12)*d(7)+c(6)*d(19)+c(19)*d(6)+c(10)*d(17)+ & + & c(17)*d(10) + s(44)=c(6)*d(13)+c(13)*d(6)+c(7)*d(17)+c(17)*d(7)+c(10)*d(19)+ & + & c(19)*d(10) + s(45)=c(8)*d(11)+c(11)*d(8)+c(5)*d(14)+c(14)*d(5) + s(46)=c(9)*d(11)+c(11)*d(9)+c(5)*d(15)+c(15)*d(5) + s(47)=c(8)*d(12)+c(12)*d(8)+c(6)*d(16)+c(16)*d(6) + s(48)=c(10)*d(12)+c(12)*d(10)+c(6)*d(17)+c(17)*d(6) + s(49)=c(10)*d(13)+c(13)*d(10)+c(7)*d(19)+c(19)*d(7) + s(50)=c(9)*d(13)+c(13)*d(9)+c(7)*d(18)+c(18)*d(7) + s(51)=c(8)*d(13)+c(13)*d(8)+c(7)*d(20)+c(20)*d(7)+c(9)*d(19)+ & + & c(19)*d(9)+c(10)*d(18)+c(18)*d(10) + s(52)=c(10)*d(11)+c(11)*d(10)+c(5)*d(20)+c(20)*d(5)+c(9)*d(14)+ & + & c(14)*d(9)+c(8)*d(15)+c(15)*d(8) + s(53)=c(9)*d(12)+c(12)*d(9)+c(6)*d(20)+c(20)*d(6)+c(10)*d(16)+ & + & c(16)*d(10)+c(8)*d(17)+c(17)*d(8) + s(54)=c(5)*d(17)+c(17)*d(5)+c(6)*d(15)+c(15)*d(6)+c(14)*d(10)+ & + & d(14)*c(10)+c(9)*d(16)+c(16)*d(9)+c(8)*d(20)+c(20)*d(8) + s(55)=c(5)*d(19)+c(19)*d(5)+c(7)*d(14)+c(14)*d(7)+c(10)*d(15)+ & + & c(15)*d(10)+c(9)*d(20)+c(20)*d(9)+c(8)*d(18)+c(18)*d(8) + s(56)=c(6)*d(18)+c(18)*d(6)+c(7)*d(16)+c(16)*d(7)+c(10)*d(20)+ & + & c(20)*d(10)+c(9)*d(17)+c(17)*d(9)+c(8)*d(19)+c(19)*d(8) + ! end of d - f + if(iff1.eq.3.or.iff2.eq.3) return + ! f-f - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1710 + s(57)=c(11)*d(11) + s(58)=c(12)*d(12) + s(59)=c(13)*d(13) + s(60)=c(11)*d(12)+c(12)*d(11)+c(14)*d(16)+c(16)*d(14) + s(61)=c(11)*d(13)+c(13)*d(11)+c(15)*d(18)+c(18)*d(15) + s(62)=c(12)*d(13)+c(13)*d(12)+c(17)*d(19)+c(19)*d(17) + s(63)=c(11)*d(14)+c(14)*d(11) + s(64)=c(11)*d(15)+c(15)*d(11) + s(65)=c(13)*d(18)+c(18)*d(13) + s(66)=c(13)*d(19)+c(19)*d(13) + s(67)=c(12)*d(17)+d(12)*c(17) + s(68)=c(12)*d(16)+c(16)*d(12) + s(69)=c(11)*d(16)+c(16)*d(11)+c(14)*d(14) + s(70)=c(11)*d(18)+c(18)*d(11)+c(15)*d(15) + s(71)=c(13)*d(15)+c(15)*d(13)+c(18)*d(18) + s(72)=c(13)*d(17)+c(17)*d(13)+c(19)*d(19) + s(73)=c(12)*d(14)+c(14)*d(12)+c(16)*d(16) + s(74)=c(12)*d(19)+c(19)*d(12)+c(17)*d(17) + s(75)=c(11)*d(17)+c(17)*d(11)+c(14)*d(20)+c(20)*d(14)+ & + & c(15)*d(16)+c(16)*d(15) + s(76)=c(11)*d(19)+c(19)*d(11)+c(20)*d(15)+c(15)*d(20)+ & + & c(14)*d(18)+c(18)*d(14) + s(77)=c(12)*d(18)+c(18)*d(12)+c(16)*d(19)+c(19)*d(16)+ & + & c(17)*d(20)+c(20)*d(17) + s(78)=c(13)*d(14)+c(14)*d(13)+c(15)*d(19)+c(19)*d(15)+ & + & c(18)*d(20)+c(20)*d(18) + s(79)=c(12)*d(15)+c(15)*d(12)+c(14)*d(17)+c(17)*d(14)+ & + & c(16)*d(20)+c(20)*d(16) + s(80)=c(13)*d(16)+c(16)*d(13)+c(17)*d(18)+c(18)*d(17)+ & + & c(19)*d(20)+c(20)*d(19) + s(81)=c(11)*d(20)+c(20)*d(11)+c(14)*d(15)+c(15)*d(14) + s(82)=c(12)*d(20)+c(20)*d(12)+c(16)*d(17)+c(17)*d(16) + s(83)=c(13)*d(20)+c(20)*d(13)+c(18)*d(19)+c(19)*d(18) + s(84)=c(14)*d(19)+c(19)*d(14)+c(15)*d(17)+c(17)*d(15)+ & + & c(16)*d(18)+c(18)*d(16)+c(20)*d(20) + ! end of f - f + return +end subroutine prod + + +!! --------------------------------------------------------------[SAW1710]- +! made pure, made explicit +pure subroutine lmnpre(l,m,n,lmnexp,lmnfak) + implicit none + integer,intent(in) :: l,m,n + real(wp), intent(out) :: lmnfak + integer,intent(out) :: lmnexp + integer :: lh,mh,nh + real(wp), parameter :: dftr(7) = & + & (/1.d0,1.d0,3.d0,15.d0,105.d0,945.d0,10395.d0/) + lmnfak=0 + lmnexp=0 + if (mod(l,2).ne.0) return + lh=l/2 + if (mod(m,2).ne.0) return + mh=m/2 + if (mod(n,2).ne.0) return + nh=n/2 + lmnexp=lh+mh+nh + lmnfak=dftr(lh+1)*dftr(mh+1)*dftr(nh+1) + return +end subroutine lmnpre + + +!! --------------------------------------------------------------[SAW1710]- +! made explicit, changed data in parameter, made elemental +elemental function olap2(l,m,n,arg,gama) + implicit none + integer,intent(in) :: l,m,n + real(wp), intent(in) :: arg,gama + real(wp) :: olap2 + integer :: lh,mh,nh + real(wp), parameter :: dftr(7) = & + & (/1.d0,1.d0,3.d0,15.d0,105.d0,945.d0,10395.d0/) + if (mod(l,2).ne.0) goto 50 + lh=l/2 + if (mod(m,2).ne.0) goto 50 + mh=m/2 + if (mod(n,2).ne.0) goto 50 + nh=n/2 + olap2=arg*gama**(lh+mh+nh)*dftr(lh+1)*dftr(mh+1)*dftr(nh+1) + return + 50 olap2=0.d0 + return +end function olap2 + +!! --------------------------------------------------------------[SAW1710]- +! removed implicit, changed data in parameter, made elemental +elemental function olap(l,m,n,gama) + implicit none + integer,intent(in) :: l,m,n + real(wp), intent(in) :: gama + real(wp) :: olap + real(wp) :: gm + integer :: lh,mh,nh + real(wp), parameter :: dftr(7) = & + & (/1.d0,1.d0,3.d0,15.d0,105.d0,945.d0,10395.d0/) + real(wp), parameter :: pi = 3.1415926535897932384626433832795029d0 + real(wp), parameter :: thehlf = 0.5d0 + gm=1.0d0/gama + if (mod(l,2).ne.0) goto 50 + lh=l/2 + if (mod(m,2).ne.0) goto 50 + mh=m/2 + if (mod(n,2).ne.0) goto 50 + nh=n/2 + olap=(dsqrt(pi*gm))**3*(thehlf*gm)**(lh+mh+nh)* & + & dftr(lh+1)*dftr(mh+1)*dftr(nh+1) + return + 50 olap=0.d0 + return +end function olap + +!! --------------------------------------------------------------[SAW1710]- +! removed implicit, made pure +pure subroutine opab1(l,m,n,ga,v,d) + ! electronic part of dipole moment + implicit none + integer,intent(in) :: l,m,n + real(wp), intent(in) :: ga,d(*) + real(wp), intent(out) :: v(*) + integer :: i + real(wp) :: g(4) + g(1)=olap(l,m,n,ga) + g(2)=olap(l+1,m,n,ga) + g(3)=olap(l,m+1,n,ga) + g(4)=olap(l,m,n+1,ga) + v(1)=g(2)+d(1)*g(1) + v(2)=g(3)+d(2)*g(1) + v(3)=g(4)+d(3)*g(1) + do i=1,3 + v(i)=-v(i) + enddo + return +end subroutine opab1 + +!! --------------------------------------------------------------[SAW1710]- +! made explicit, made pure +pure subroutine opab4(l,m,n,ga,v,d) + ! electronic part of second moment + implicit none + integer,intent(in) :: l,m,n + real(wp), intent(in) :: ga,d(*) + real(wp), intent(out) :: v(*) + real(wp) :: g(10) + integer :: i + g(1)=olap(l,m,n,ga) + g(2)=olap(l+1,m,n,ga) + g(3)=olap(l,m+1,n,ga) + g(4)=olap(l,m,n+1,ga) + g(5)=olap(l+2,m,n,ga) + g(6)=olap(l,m+2,n,ga) + g(7)=olap(l,m,n+2,ga) + g(8)=olap(l+1,m+1,n,ga) + g(9)=olap(l+1,m,n+1,ga) + g(10)=olap(l,m+1,n+1,ga) + v(1)=g(5)+2.d0*d(1)*g(2)+d(1)**2*g(1) + v(2)=g(6)+2.d0*d(2)*g(3)+d(2)**2*g(1) + v(3)=g(7)+2.d0*d(3)*g(4)+d(3)**2*g(1) + v(4)=g(8)+d(1)*g(3)+d(2)*g(2)+d(1)*d(2)*g(1) + v(5)=g(9)+d(1)*g(4)+d(3)*g(2)+d(1)*d(3)*g(1) + v(6)=g(10)+d(2)*g(4)+d(3)*g(3)+d(2)*d(3)*g(1) + do i=1,6 + v(i)=-v(i) + enddo + return +end subroutine opab4 + +!! --------------------------------------------------------------[SAW1710]- +! made explict, made pure +pure subroutine opac3(l,m,n,ga,v,d) + ! electronic part of charge density + implicit none + integer,intent(in) :: l,m,n + real(wp), intent(in) :: ga,d(3) + real(wp), intent(out) :: v(1) + real(wp) :: dd(3),ddsq + integer :: i + ddsq=0.d0 + do i=1,3 + dd(i)=-d(i) + ddsq=ddsq+d(i)**2 + enddo + v(1) =1.d0 + if(l.ne.0) v(1)=dd(1)**l + if(m.ne.0) v(1)=v(1)*dd(2)**m + if(n.ne.0) v(1)=v(1)*dd(3)**n + v(1)=v(1)*dexp(-ga*ddsq) + return +end subroutine opac3 + +!! --------------------------------------------------------------[SAW1710]- +! made explicit, made pure +pure subroutine opad1(l,m,n,ga,v,d) + ! electronic part of overlap + implicit none + integer,intent(in) :: l,m,n + real(wp), intent(in) :: ga,d(3) + real(wp), intent(out) :: v(1) + v(1)=olap(l,m,n,ga) + return +end subroutine opad1 + + +!! --------------------------------------------------------------[SAW1710]- +! changed do loops +subroutine opap4(l,m,n,gc,v,rc) + ! + ! this subroutine calculates the expectation-values of + ! -p**2/2 (kinetic energy) and p**4/(8*cl**2) (relativistic + ! correction of the kinetic energy) + ! + ! the program works correct only for s and p and the irreducible + ! linear combinations of the primitive d,f,... functions + ! + ! written by s. brode in april,1984 + ! + implicit real(wp) (a-h,o-z) + logical :: lodd,modd,nodd,leven,meven,neven + real(wp) :: v(1),rc(3),rca(3),rcb(3),srcab(3),o(25) + ! ----- ra,rb centres ga,gb exponents ia,ib monom indices + common /abfunc/ ra(3),rb(3),ga,gb,ia,ib + real(wp), parameter :: cl = 137.03604d0 + real(wp), parameter :: dftr(7) = & + & (/1.d0,1.d0,3.d0,15.d0,105.d0,945.d0,10395.d0/) + real(wp), parameter :: pi=3.1415926535897932384626433832795029d0 + real(wp), parameter :: a0=0.0d0,a1=1.0d0,a2=2.0d0,a4=4.0d0,a8=8.0d0 + integer,parameter :: lmni(5) = (/0,3*1,6*2,10*3,15*4/) + ! + ! define the overlap-integral (in line function definition) + ! + !!!ola(ll,mm,nn)=gch**(ll/2+mm/2+nn/2)*dftr(ll/2+1)*dftr(mm/2+1)* & + !!! & dftr(nn/2+1) + ! + ! some previous calculations + ! + ! do 10 i=1,25 + ! 10 o(i)=a0 + o(1:7)=0 + gch=a1/(a2*gc) + ropigc=dsqrt(pi/gc)*pi/gc + lodd=mod(l,2).eq.1 + leven=.not.lodd + modd=mod(m,2).eq.1 + meven=.not.modd + nodd=mod(n,2).eq.1 + neven=.not.nodd + ga2=ga*a2 + gb2=gb*a2 + gab4=ga2*gb2 + prcaa=a0 + prcbb=a0 + do i=1,3 + rca(i)=rc(i)-ra(i) + rcb(i)=rc(i)-rb(i) + srcab(i)=rca(i)+rcb(i) + prcaa=prcaa+rca(i)*rca(i) + prcbb=prcbb+rcb(i)*rcb(i) + enddo + zwa=ga2*prcaa-dble(2*lmni(ia)+3) + zwb=gb2*prcbb-dble(2*lmni(ib)+3) + ! + ! calculation of the overlap-integrals + ! + if(lodd.or.modd.or.nodd) go to 110 + o( 1)=ola(l ,m ,n ) + o( 5)=ola(l+2,m ,n ) + o( 6)=ola(l ,m+2,n ) + o( 7)=ola(l ,m ,n+2) + ! o(11)=ola(l+4,m ,n ) + ! o(12)=ola(l ,m+4,n ) + ! o(13)=ola(l ,m ,n+4) + ! o(20)=ola(l+2,m+2,n ) + ! o(21)=ola(l+2,m ,n+2) + ! o(22)=ola(l ,m+2,n+2) + 110 continue + ! + if(leven.or.modd.or.nodd) go to 120 + o( 2)=ola(l+1,m ,n ) + ! o( 8)=ola(l+3,m ,n ) + ! o(14)=ola(l+1,m+2,n ) + ! o(15)=ola(l+1,m ,n+2) + 120 continue + ! + if(lodd.or.meven.or.nodd) go to 130 + o( 3)=ola(l ,m+1,n ) + ! o( 9)=ola(l ,m+3,n ) + ! o(16)=ola(l+2,m+1,n ) + ! o(17)=ola(l ,m+1,n+2) + 130 continue + ! + if(lodd.or.modd.or.neven) go to 140 + o( 4)=ola(l ,m ,n+1) + ! o(10)=ola(l ,m ,n+3) + ! o(18)=ola(l+2,m ,n+1) + ! o(19)=ola(l ,m+2,n+1) + 140 continue + ! + ! if(leven.or.meven.or.nodd) go to 150 + ! o(23)=ola(l+1,m+1,n ) + ! 150 continue + ! if(leven.or.modd.or.neven) go to 160 + ! o(24)=ola(l+1,m ,n+1) + ! 160 continue + ! if(lodd.or.meven.or.neven) go to 170 + ! o(25)=ola(l ,m+1,n+1) + ! 170 continue + ! + ! calculation of kinetic energy + ! + p2=-ga*(zwa*o(1) & + & +ga2*(a2*(rca(1)*o(2)+rca(2)*o(3)+rca(3)*o(4)) & + & +o(5)+o(6)+o(7))) + + + ! calculation of relativistic correction + + ! p4=(gab4*(o(11)+o(12)+o(13)+a2*(o(20)+o(21)+o(22)) + ! . +a2*(srcab(1)*(o( 8)+o(16)+o(18)) + ! . +srcab(2)*(o(14)+o( 9)+o(19)) + ! . +srcab(3)*(o(15)+o(17)+o(10))) + ! . +a4*(rca(1)*rcb(1)*o(5) + ! . +rca(2)*rcb(2)*o(6) + ! . +rca(3)*rcb(3)*o(7) + ! . +a2*((rca(1)*rcb(2)+rca(2)+rcb(1))*o(23) + ! . +(rca(1)*rcb(3)+rca(3)+rcb(1))*o(24) + ! . +(rca(2)*rcb(3)+rca(3)+rcb(2))*o(25)))) + ! . +(gb2*zwa+ga2*zwb)*(o(5)+o(6)+o(7)) + ! . +a2*(gb2*zwa*(rcb(1)*o(2)+rcb(2)*o(3)+rcb(3)*o(4)) + ! . +ga2*zwb*(rca(1)*o(2)+rca(2)*o(3)+rca(3)*o(4))) + ! . +zwa*zwb*o(1)) + ! . *(-gab4)/(a8*cl**2) + + ! store p2 and p4 to the correct storage-locations + + v(1)=p2*ropigc + ! v(2)=p4*ropigc + ! v(3)=v(1)+v(2) + return + +contains + + real(wp) function ola(ll,mm,nn) + integer, intent(in) :: ll,mm,nn + ola=gch**(ll/2+mm/2+nn/2)*dftr(ll/2+1)*dftr(mm/2+1)*dftr(nn/2+1) + end function ola + +end subroutine opap4 + + +!subroutine opaa0(l,m,n,ga,v,d) +! ! electronic part of potential +! implicit real(wp)(a-h,o-z) +! dimension v(*),d(*),fnu(7),dp(3) +! real(wp), parameter :: fn(7) = & +! & (/1.d0,1.d0,2.d0,6.d0,24.d0,120.d0,720.d0/) +! real(wp), parameter :: fd(7) = (/1.d0,1.d0,0.5d0,0.1666666666667d0, & +! & 0.416666666667d-1,.833333333333d-2,.1388888888889d-2/) +! real(wp), parameter :: pi = 3.1415926535897932384626433832795029d0 +! itot=l+m+n +! itoth=itot/2 +! pre=(2.d0*pi/ga)*fn(l+1)*fn(m+1)*fn(n+1) +! !!!if(2*itoth-itot) 1,2,1 +! ival=2*itoth-itot +! if(ival.eq.0) then +! goto 2 +! else +! goto 1 +! endif +! 1 pre=-pre +! 2 del=.25d0/ga +! x=ga*(d(1)**2+d(2)**2+d(3)**2) +! xx=2.d0*x +! call fmc(6,x,expmx,fmch) +! fnu(7)=fmch +! fnu(6)=(expmx+xx*fnu(7))/11.d0 +! fnu(5)=(expmx+xx*fnu(6))/9.d0 +! fnu(4)=(expmx+xx*fnu(5))/7.d0 +! fnu(3)=(expmx+xx*fnu(4))/5.d0 +! fnu(2)=(expmx+xx*fnu(3))/3.d0 +! fnu(1)=expmx+xx*fnu(2) +! dp(1)=1.d0 +! dp(2)=del +! dp(3)=del**2 +! v(1)= - pre*aainer(0,0,0,l,m,n,d,dp,fnu,fn,fd) +! return +!end subroutine opaa0 + + +!double precision function aainer(r,s,t,l,m,n,d,dp,fnu,fn,fd) +! implicit real(wp)(a-h,o-z) +! dimension d(*),dp(*),fnu(*),fn(*),fd(*) +! integer r,s,t,u,v,w,uvwt,rstt,uvwth +! rstt=r+s+t +! lmnt=l+m+n +! lmnrst=rstt+lmnt +! lh=l/2 +! mh=m/2 +! nh=n/2 +! aainer=0.0d0 +! last1 = lh+1 +! do ii=1,last1 +! i = ii-1 +! il=l-2*i+1 +! ilr=r+il +! ilrh=(ilr-1)/2 +! if (ilr.gt.9 .or. i+1.gt.9 .or. il.gt.9) stop 'aainer: ilr' +! fi=fn(ilr)*fd(il)*fd(i+1) +! last2 = mh+1 +! do jj=1,last2 +! j = jj-1 +! jm=m-2*j+1 +! jms=s+jm +! jmsh=(jms-1)/2 +! if (jms.gt.9 .or. j+1.gt.9 .or. jm.gt.9) stop 'aainer: ilr' +! fij=fn(jms)*fd(jm)*fd(j+1)*fi +! last3 = nh+1 +! do kk=1,last3 +! k = kk-1 +! kn=n-2*k+1 +! knt=t+kn +! knth=(knt-1)/2 +! ijkt=i+j+k +! if (knt.gt.9 .or. k+1.gt.9 .or. kn.gt.9) stop 'aainer: ilr' +! fijk=fn(knt)*fd(kn)*fd(k+1)*dp(ijkt+1)*fij +! lmrsij=lmnrst-2*ijkt +! last4 = ilrh+1 +! do iu=1,last4 +! u = iu-1 +! ilru=ilr-2*u +! if (u+1.gt.9 .or. ilru.gt.9) stop 'aainer: ilru' +! fu=fd(u+1)*fd(ilru) +! !!!if(dabs(d(1))-0.0000001d0)10,10,11 +! dval=dabs(d(1))-0.0000001d0 +! if(dval.le.0d0) then +! goto 10 +! else +! goto 11 +! endif +! !!!10 if(ilru-1) 12,12, 3 +! 10 ival=ilru-1 +! if(ival.le.0) then +! goto 12 +! else +! goto 3 +! endif +! 11 fu=fu*d(1)**(ilru-1) +! 12 last5 = jmsh+1 +! do iv=1,last5 +! v = iv-1 +! jmsv=jms-2*v +! fuv=fu*fd(v+1)*fd(jmsv) +! !!!if(dabs(d(2))-0.0000001d0)20,20,21 +! dval=dabs(d(2))-0.0000001d0 +! if(dval.le.0d0) then +! goto 20 +! else +! goto 21 +! endif +! !!!20 if(jmsv-1) 22,22,2 +! 20 ival=jmsv-1 +! if(ival.le.0) then +! goto 22 +! else +! goto 2 +! endif +! 21 fuv=fuv*d(2)**(jmsv-1) +! 22 last6 = knth+1 +! do iw=1,last6 +! w = iw-1 +! kntw=knt-2*w +! fuvw=fuv*fd(w+1) *fd(kntw) +! !!!if(dabs(d(3))-0.0000001d0)30,30,31 +! dval=dabs(d(3))-0.0000001d0 +! if(dval.le.0d0) then +! goto 30 +! else +! goto 31 +! endif +! !!!30 if(kntw-1) 32,32,1 +! 30 ival=kntw-1 +! if(ival.le.0) then +! goto 32 +! else +! goto 1 +! endif +! 31 fuvw=fuvw*d(3)**(kntw-1) +! 32 uvwt=u+v+w +! uvwth=uvwt/2 +! !!!if(2*uvwth-uvwt) 33,34,33 +! ival=2*uvwth-uvwt +! if(ival.eq.0) then +! goto 34 +! else +! goto 33 +! endif +! 33 fuvw=-fuvw +! 34 nuindx=lmrsij-uvwt +! fuvw=fuvw*fnu(nuindx +1)*dp(uvwt+1) +! aainer=fijk*fuvw+aainer +! 1 continue +! enddo +! 2 continue +! enddo +! 3 continue +! enddo +! enddo +! enddo +! enddo +! +! return +!end function aainer + + +subroutine fmc(mvar,xvar,expmx,fmch) + implicit real(wp) (a-h,o-z) + m=mvar + x=xvar + !mk ---- prevent underflow error messages + expmx=0.d0 + if(x.lt.50.d0) expmx=dexp(-x) + ! expmx=dexp(-x) + ! + !!!if(x -20.d0) 10,10,20 + dval=x-20d0 + if(dval.le.0d0) then + goto 10 + else + goto 20 + endif + ! x less than or equal to 20.0. + 10 a=m + a=a+0.5d0 + term=1.d0/a + ptlsum=term + do i=2,60 + a=a+1.d0 + term=term*x/a + ptlsum=ptlsum+term + !!!if (term/ptlsum- 1.d-8 ) 12, 11, 11 + dval=term/ptlsum-1.d-8 + if(dval.lt.0d0) then + goto 12 + else + goto 11 + endif + 11 continue + enddo + write(6,'(1x,a,i6,d16.9)') 'no convergence for fmch',m,x + 12 fmch=.5d0*ptlsum*expmx + return + ! x greater than 20.0. + 20 a=m + xd=1.d0/x + b=a+0.5d0 + a=a-0.5d0 + approx=.886226925428d0*(dsqrt(xd)*xd**m) + !!!if (m) 21, 23, 21 + if(m.ne.0d0) then + do i=1,m + b=b-1.d0 + approx=approx*b + enddo + endif + fimult=.5d0*expmx*xd + fiprop=fimult/approx + term=1.d0 + ptlsum=term + notrms=int(x) + notrms=notrms+m + do i=2,notrms + term=term*a*xd + ptlsum=ptlsum+term + !!!if (dabs(term*fiprop/ptlsum) - 1.d-8 ) 25,25,24 + dval=dabs(term*fiprop/ptlsum)-1.d-8 + if(dval.le.0d0) then + goto 25 + else + goto 24 + endif + 24 continue + a=a-1.d0 + enddo + + write(6,'(1x,a,i6,d16.9)') 'no convergence for fmch',m,x + 25 fmch=approx-fimult*ptlsum + return +end subroutine fmc + + subroutine opaa0(l,m,n,ga,v,d) +! electronic part of potential + implicit real*8(a-h,o-z) + dimension v(*),d(*),fnu(7),fn(7),fd(7),dp(3) + data fn /1.d0,1.d0,2.d0,6.d0,24.d0,120.d0, & + & 720.d0/ + data fd/1.d0,1.d0,0.5d0,0.1666666666667d0,0.416666666667d-1, & + & .833333333333d-2,.1388888888889d-2/ + data pi/3.1415926535897932384626433832795029d0/ + itot=l+m+n + itoth=itot/2 + pre=(2.d0*pi/ga)*fn(l+1)*fn(m+1)*fn(n+1) + if (modulo(itot, 2) == 1) pre = -pre + del=.25d0/ga + x=ga*(d(1)**2+d(2)**2+d(3)**2) + xx=2.d0*x + call fmc (6,x,expmx,fmch) ! 50 % of CPU + fnu(7)=fmch + fnu(6)=(expmx+xx*fnu(7))/11.d0 + fnu(5)=(expmx+xx*fnu(6))/9.d0 + fnu(4)=(expmx+xx*fnu(5))/7.d0 + fnu(3)=(expmx+xx*fnu(4))/5.d0 + fnu(2)=(expmx+xx*fnu(3))/3.d0 + fnu(1)=expmx+xx*fnu(2) + dp(1)=1.d0 + dp(2)=del + dp(3)=del**2 + v(1)= pre*aainer(l,m,n,d,dp,fnu,fn,fd) + return + end + + + double precision function aainer(l,m,n,d,dp,fnu,fn,fd) + implicit real*8(a-h,o-z) + dimension d(*),dp(*),fnu(*),fn(*),fd(*) + integer u,v,w,uvwt,rstt,uvwth +! rstt=r+s+t + lmnt=l+m+n +! lmnrst=rstt+lmnt + lmnrst= lmnt + lh=l/2 + mh=m/2 + nh=n/2 + aainer=0.0d0 + last1 = lh+1 + do ii=1,last1 + i = ii-1 + il=l-2*i+1 +! ilr=r+il + ilr=il + ilrh=(ilr-1)/2 +!sg if (ilr.gt.9 .or. i+1.gt.9 .or. il.gt.9) stop 'aainer: ilr' + fi=fn(ilr)*fd(il)*fd(i+1) + last2 = mh+1 + do jj=1,last2 + j = jj-1 + jm=m-2*j+1 +! jms=s+jm + jms=jm + jmsh=(jms-1)/2 +!sg if (jms.gt.9 .or. j+1.gt.9 .or. jm.gt.9) stop 'aainer: ilr' + fij=fn(jms)*fd(jm)*fd(j+1)*fi + last3 = nh+1 + do kk=1,last3 + k = kk-1 + kn=n-2*k+1 +! knt=t+kn + knt=kn + knth=(knt-1)/2 + ijkt=i+j+k +!sg if (knt.gt.9 .or. k+1.gt.9 .or. kn.gt.9) stop 'aainer: ilr' + fijk=fn(knt)*fd(kn)*fd(k+1)*dp(ijkt+1)*fij + lmrsij=lmnrst-2*ijkt + last4 = ilrh+1 + do iu=1,last4 + u = iu-1 + ilru=ilr-2*u +!sg if (u+1.gt.9 .or. ilru.gt.9) stop 'aainer: ilru' + fu=fd(u+1)*fd(ilru) + if(dabs(d(1))-0.0000001d0.le.0.0) then + if(ilru-1.gt.0) cycle + else + fu=fu*d(1)**(ilru-1) + end if + last5 = jmsh+1 + do iv=1,last5 + v = iv-1 + jmsv=jms-2*v + fuv=fu*fd(v+1)*fd(jmsv) + if(dabs(d(2))-0.0000001d0.le.0d0) then + if(jmsv-1.gt.0) cycle + else + fuv=fuv*d(2)**(jmsv-1) + end if + last6 = knth+1 + do iw=1,last6 + w = iw-1 + kntw=knt-2*w + fuvw=fuv*fd(w+1) *fd(kntw) + if(dabs(d(3))-0.0000001d0.le.0d0) then + if(kntw-1.gt.0) cycle + else + fuvw=fuvw*d(3)**(kntw-1) + end if + uvwt=u+v+w + uvwth=uvwt/2 + if (modulo(uvwt, 2) == 1) fuvw = -fuvw + nuindx=lmrsij-uvwt + fuvw=fuvw*fnu(nuindx+1)*dp(uvwt+1) + aainer=fijk*fuvw+aainer + enddo + enddo + enddo + enddo + enddo + enddo + return + end + + +end module xtb_intpack diff --git a/src/qmhelpers/local.f90 b/src/qmhelpers/local.f90 new file mode 100644 index 00000000..d4d8cce9 --- /dev/null +++ b/src/qmhelpers/local.f90 @@ -0,0 +1,1248 @@ +! Copyright (C) 2017-2020 Stefan Grimme, Sebastian Ehlert (xtb) +! Copyright (C) 2024 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 +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! crest is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with crest. If not, see . +! +! Routines were adapted from the xtb code (github.com/grimme-lab/xtb) +! under the Open-source software LGPL-3.0 Licencse. +!================================================================================! +module mo_localize + use crest_parameters + use lopt_mod,only:lopt + use strucrd,only:i2e + use wiberg_mayer + implicit none + private + + public :: local + + logical,parameter :: pr_local = .true. + +!========================================================================================! +!========================================================================================! +contains !> MODULE PROCEDURES START HERE +!========================================================================================! +!========================================================================================! + + subroutine local(nat,at,nbf,nao,ihomo,xyz,z, & + & focc,s,p,cmo,eig, & + & aoat2,aoat,nprim,alp,lao,cont, & + & nprot,protxyz ) + implicit none + !> INPUT + integer,intent(in) :: nat !> number of atoms + integer,intent(in) :: at(nat) !> atomic numbers for each atom + integer,intent(in) :: nbf !> number of basis functions + integer,intent(in) :: nao !> number of AOs + integer,intent(in) :: ihomo !> number of highest occ MO + real(wp),intent(in) :: xyz(3,nat) !> molecular coordinates (Bohr) + real(wp),intent(in) :: z(nat) !> nuclear charge + real(wp),intent(in) :: focc(nao) !> occupations + real(wp),intent(in) :: s(nao,nao) !> overlap S + real(wp),intent(in) :: p(nao,nao) !> density matrix P + real(wp),intent(in) :: cmo(nao,nao) !> MO coefficients C + real(wp),intent(in) :: eig(nao) !> orbital energies + !real(wp),intent(in) :: q(nat) !> atomic charges + + integer,intent(inout) :: aoat2(nao) !> mapping AO -> atom number + integer,intent(in) :: aoat(nbf) !> mapping BF -> atom + integer,intent(in) :: nprim(nbf) !> mappinf BF -> primitive number + real(wp),intent(in) :: alp(9*nbf) !> mapping primitive -> primitive exponent + integer,intent(in) :: lao(nbf) !> mapping BF -> azimudal quantum number l + real(wp),intent(in) :: cont(9*nbf) !> mapping primitive -> contraction coeffient + + !> OUTPUT + integer,intent(out) :: nprot + real(wp),intent(out),allocatable :: protxyz(:,:) + + !> LOCAL + real(wp),allocatable :: op(:,:) + real(wp),allocatable :: oc(:,:,:) + real(wp),allocatable :: cca(:) + real(wp),allocatable :: dip2(:) + real(wp),allocatable :: qua(:,:) + real(wp),allocatable :: dip(:,:) + real(wp),allocatable :: ecent(:,:),qcent(:,:),ecent2(:,:) + real(wp),allocatable :: d(:,:) + real(wp),allocatable :: f(:) + real(wp),allocatable :: eiga(:) + real(wp),allocatable :: xcen(:) + real(wp),allocatable :: qmo(:,:) + real(wp),allocatable :: tmpq(:,:) + real(wp),allocatable :: rr(:) + real(wp),allocatable :: wbo(:,:) + real(wp),allocatable :: xyztmp(:,:) + real(sp),allocatable :: rklmo(:,:) + integer,allocatable :: ind(:) + integer,allocatable :: lneigh(:,:) + integer,allocatable :: aneigh(:,:) + + integer :: mo,n,i,j,k,ii,jj,imem(nat),idum,isc,iso,nop,pilist(ihomo) + integer :: lamdu1,lamdu2,imo1,imo2,iso1,iso2,ij,jdum,maxlp,maxpi,is1 + integer :: ilumo,klev,nlev,nn,m,ldum(ihomo),sigrel(3,ihomo),npi,is2 + integer :: i1,i2,i3,is3,ipi,piset(nat),ncyc,j1,j2,smo,pmo,nl,nm,nr + integer :: imo,new,nao2 + real(wp) :: dd,dum,det,pp,dtot(3),diptot,thr,t0,w0,t1,w1,vec1(3),v(6) + real(wp) :: enlumo,enhomo,qhl(nat,2),ps,efh,efl,r1,r2,pithr,vec2(3) + real(wp) :: praxis(3,3),aa,bb,cc + character(len=80) :: atmp + character(len=5) :: lmostring(4) + data lmostring/'sigma','LP ','pi ','delpi'/ + !data lmostring/'σ ','LP ','π ','δπ '/ + logical l1,l2,l3,flip + integer LWORK,IWORK,LIWORK,INFO + integer :: iscreen,icoord,ilmoi,icent ! file handles + + !> BLAS + external dgemm + +! if (pr_local) then +! write (*,*) +! write (*,*) 'localization/xTB-IFF output generation' +! end if + n = ihomo + ilumo = n+1 + + +! if (ilumo .gt. nao) then +! enlumo = 1000.d0 +! else +! enlumo = eig(ilumo) +! end if +! efh = eig(n) +! efl = enlumo +! +! !> HOMO/LUMO populations for xTB FF (CT terms) +! qhl = 0 +! thr = 0.01 +! +! if (ihomo .eq. 0) then +! !> the HOMO is non-existing, place at unrealistically low energy +! enhomo = -9999.999999990d0 +! qhl(1:nat,1) = 0.0d0 +! else +! klev = 0 +! enhomo = 0 +! do nlev = ihomo,1,-1 ! loop over highest levels +! if (efh-eig(nlev) .gt. thr) exit +! klev = klev+1 +! enhomo = enhomo+eig(nlev) +! do i = 1,nao +! ii = aoat2(i) +! do j = 1,i-1 +! jj = aoat2(j) +! ps = s(j,i)*cmo(j,nlev)*cmo(i,nlev) +! qhl(ii,1) = qhl(ii,1)+ps +! qhl(jj,1) = qhl(jj,1)+ps +! end do +! ps = s(i,i)*cmo(i,nlev)*cmo(i,nlev) +! qhl(ii,1) = qhl(ii,1)+ps +! end do +! end do +! dum = 1.0_wp/float(klev) +! enhomo = enhomo*dum +! qhl(1:nat,1) = qhl(1:nat,1)*dum +! if (pr_local) write (*,*) 'averaging CT terms over ',klev,' occ. levels' +! end if +! +! klev = 0 +! enlumo = 0 +! do nlev = ilumo,nao !> loop over highest levels +! if (eig(nlev)-efl .gt. thr) exit +! klev = klev+1 +! enlumo = enlumo+eig(nlev) +! do i = 1,nao +! ii = aoat2(i) +! do j = 1,i-1 +! jj = aoat2(j) +! ps = s(j,i)*cmo(j,nlev)*cmo(i,nlev) +! qhl(ii,2) = qhl(ii,2)+ps +! qhl(jj,2) = qhl(jj,2)+ps +! end do +! ps = s(i,i)*cmo(i,nlev)*cmo(i,nlev) +! qhl(ii,2) = qhl(ii,2)+ps +! end do +! end do +! dum = 1.0d0/float(klev) +! enlumo = enlumo*dum +! qhl(1:nat,2) = qhl(1:nat,2)*dum +! if (pr_local) write (*,*) 'averaging CT terms over ',klev,' virt. levels' +! if (ilumo .gt. nao) then +! enlumo = 1000.0d0 +! qhl(1:nat,2) = 0 +! end if + + allocate (cca(nao*nao),xcen(n),lneigh(4,n),aneigh(2,n)) + allocate (d(n,n),ecent(n,4),eiga(n),qcent(n,3),ecent2(n,4)) + + !> do only occ. ones + cca = 0 + k = 0 + do i = 1,n + k = k+1 + eiga(k) = eig(i) + do j = 1,nao + cca(j+(k-1)*nao) = cmo(j,i) + end do + end do + + nop = 3 +!>--- dipole integrals + allocate (dip2(nao*nao),dip(nbf*(nbf+1)/2,3),op(n*(n+1)/2,nop)) + call Dints(nat,nbf,xyz,dip(:,1),dip(:,2),dip(:,3), & + & aoat,nprim,alp,lao,cont) + + nao2 = nao + call cao2saop(nbf,nao2,dip(:,1),lao) + call cao2saop(nbf,nao2,dip(:,2),lao) + call cao2saop(nbf,nao2,dip(:,3),lao) + do mo = 1,3 + call onetri(1,dip(:,mo),dip2,cca,nao2,n) + k = 0 + do i = 1,n + do j = 1,i + k = k+1 + op(k,mo) = dip2(j+(i-1)*n) + end do + end do + end do + +!>--- compute dipole moment core part + dtot = 0 + do i = 1,nat + dtot(1) = dtot(1)+xyz(1,i)*z(i) + dtot(2) = dtot(2)+xyz(2,i)*z(i) + dtot(3) = dtot(3)+xyz(3,i)*z(i) + end do + k = 0 + do i = 1,n + k = k+i + dtot(1) = dtot(1)+op(k,1)*focc(i) + dtot(2) = dtot(2)+op(k,2)*focc(i) + dtot(3) = dtot(3)+op(k,3)*focc(i) + end do + diptot = sqrt(dtot(1)**2+dtot(2)**2+dtot(3)**2) + if (pr_local) then + write (*,*) 'dipole moment from electron density (au)' + write (*,*) ' X Y Z ' + write (*,'(3f9.4,'' total (Debye): '',f8.3)') & + & dtot(1),dtot(2),dtot(3),diptot*2.5418_wp + end if + +!>--- do optimization (i.e. the actual localization) + if (pr_local) write (*,*) 'doing rotations ...' + call lopt(.true.,n,3,1.d-6,op,d) + + if (pr_local) write (*,*) 'doing transformations ...' +!>--- MO charge center + k = 0 + do i = 1,n + k = k+i + ecent(i,1) = -op(k,1) + ecent(i,2) = -op(k,2) + ecent(i,3) = -op(k,3) + end do + +!>--- LMO fock matrix + allocate (f(n)) + do i = 1,n + dum = 0.0d0 + do k = 1,n + dum = dum+eiga(k)*d(k,i)*d(k,i) + end do + f(i) = dum + end do + call lmosort2(n,f,d,ecent) + +!>--- get the lmos + CALL dgemm('N','N',nao2,n,n,1.D0,cca,nao2,d,n,0.D0,cmo,nao2) ! non-std BLAS + + pithr = 2.20d0 ! below is pi, large is pi deloc (typ=4) + +!>--- number of centers for each mo + allocate (qmo(nat,n)) + call mocent(nat,nao2,n,cmo,s,qmo,xcen,aoat2) + + allocate (rklmo(5,2*n)) + !if (pr_local) write (*,*) 'lmo centers(Z=2) and atoms on file ' + if (pr_local) write (*,*) 'LMO Fii/eV ncent charge center contributions...' +! if (pr_local) open (newunit=iscreen,file='xtbscreen.xyz') + allocate (tmpq(nat,n)) + tmpq(1:nat,1:n) = qmo(1:nat,1:n) + maxlp = 0 + maxpi = 0 + do i = 1,n + do j = 1,nat + imem(j) = j + end do + call lmosort(nat,n,i,imem,qmo) + idum = 1 + do j = 1,nat + if (qmo(j,i) .gt. 0.07) idum = j + end do + if (nat .eq. 1) then + jdum = 2 + else + call lmotype(nat,at,xyz,ecent(i,1),ecent(i,2),ecent(i,3), & + & imem(1),imem(2),xcen(i),.false.,pithr,jdum) + end if + rklmo(1:3,i) = ecent(i,1:3) + rklmo(4,i) = ecent(i,4) + rklmo(5,i) = real(jdum) + if (jdum .eq. 2) maxlp = i + if (jdum .eq. 3) maxpi = i + end do + qmo(1:nat,1:n) = tmpq(1:nat,1:n) + deallocate (tmpq) + + do i = 1,n + do j = 1,nat + imem(j) = j + end do + call lmosort(nat,n,i,imem,qmo) + idum = 1 + do j = 1,nat + if (qmo(j,i) .gt. 0.07) idum = j + end do + if (nat .eq. 1) then + jdum = 1 + else + call lmotype(nat,at,xyz,ecent(i,1),ecent(i,2),ecent(i,3), & + & imem(1),imem(2),xcen(i),.true.,pithr,jdum) + end if + if (pr_local) then + write (*,'(i5,1x,a5,1x,2f7.2,3f10.5,12(i5,2x,a2,'':'',f6.2))') & + & i,lmostring(jdum),autoev*f(i),xcen(i),ecent(i,1:3), & + & (imem(j),i2e(at(imem(j))),qmo(j,i),j=1,idum) + end if + +!>--- write + LP/pi as H for protonation search + if (pr_local) then + if (jdum .gt. 1) then + if (i .eq. maxlp.or.(i .eq. maxpi.and.maxlp .eq. 0)) then +! open (newunit=icoord,file='coordprot.0') +! write (icoord,'(''$coord'')') +! do ii = 1,nat +! write (icoord,'(3F24.10,5x,a2)') xyz(1:3,ii),i2e(at(ii)) +! end do +! write (icoord,'(3F24.10,5x,a2)') ecent(i,1:3),i2e(1) +! write (icoord,'(''$end'')') +! close (icoord) +! else +! write (iscreen,*) nat+1 +! write (iscreen,*) +! do ii = 1,nat +! write (iscreen,'(a2,3F24.10)') i2e(at(ii)),xyz(1:3,ii)*autoaa +! end do +! write (iscreen,'(a2,3F24.10)') i2e(1),ecent(i,1:3)*autoaa + end if + end if + end if + + end do + + !ccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + new = 0 + if (maxval(rklmo(5,1:n)) .ge. 3.0d0) then + if (pr_local) write (*,*) 'starting deloc pi regularization ...' + + call atomneigh(n,nat,xyz,ecent,aneigh) ! the nearest atom neighbors for each LMO + call lmoneigh(n,rklmo,ecent,aneigh,lneigh) ! the nearest LMO neighbor but not LP + + ldum = -1 + do i = 1,n + if (int(rklmo(5,i)) .ge. 3.and.ldum(i) .ne. 1) then + flip = .false. + !> check for same alignment of pi orbitals + call piorient(ecent(i,:),ecent(lneigh(1,i),:),flip) + if (flip) then + ldum(i) = 1 + !> i itself is the pi DOWN counterpart and lneigh(1,i) contains the pi UP part + ldum(lneigh(1,i)) = 0 + else + ldum(i) = 0 + ldum(lneigh(1,i)) = 1 !> lneigh(1,i) contains the pi DOWN counterpart + end if + end if + end do + k = 0 + piset = 0 + do i = 1,n + if (ldum(i) .eq. 0) then + k = k+1 + pilist(k) = i !> UP plane unique pi list + piset(aneigh(1,i)) = 1 + piset(aneigh(2,i)) = 1 + end if + end do + npi = k + +!>--- determine for each deloc pi bond sigma bonds which are in the pi system + j = 0 + ldum = 0 + nn = 0 + sigrel = 0 + do i = 1,npi + ipi = pilist(i) + call irand3(i1,i2,i3) !> change order of lookup + is1 = lneigh(i1+1,ipi) + is2 = lneigh(i2+1,ipi) + is3 = lneigh(i3+1,ipi) + !> do bonds on ipi belong to the pi system? + l1 = bndcheck(nat,piset,aneigh(1,is1),aneigh(2,is1)) + l2 = bndcheck(nat,piset,aneigh(1,is2),aneigh(2,is2)) + l3 = bndcheck(nat,piset,aneigh(1,is3),aneigh(2,is3)) + if (l1.and.ldum(is1) .eq. 0) then + ldum(is1) = 1 + sigrel(1,i) = is1 + end if + if (l2.and.ldum(is2) .eq. 0) then + ldum(is2) = 1 + sigrel(2,i) = is2 + end if + if (l3.and.ldum(is3) .eq. 0) then + ldum(is3) = 1 + sigrel(3,i) = is3 + end if + end do + + ldum = 0 + do i = 1,npi + ipi = pilist(i) + do j = 1,3 + is1 = sigrel(j,i) + if (is1 .gt. 0) ldum(is1) = ldum(is1)+1 + end do + end do + + if (pr_local) write (*,*) 'thr ',pithr,'# pi deloc LMO',npi + + allocate (wbo(nat,nat)) + wbo = 0.0d0 + call get_wbo_rhf(nat,nao2,P,S,aoat2,wbo) + +!>--- now create new LMO + k = 0 + m = 0 + do i = 1,npi + pmo = pilist(i) + i1 = aneigh(1,pmo) + i2 = aneigh(2,pmo) + dd = wbo(i2,i1) + do j = 1,3 + smo = sigrel(j,i) + if (rklmo(5,smo) .ne. 1) cycle + j1 = aneigh(1,smo) + j2 = aneigh(2,smo) + dum = wbo(j2,j1) + !> take difference of pi-bond order (i.e., WBO > 1) relative to its sum (i.e., number of pi-electrons distributed among the 3 centers) + pp = abs((dum-dd)/(dum+dd-2.0d0)) + if (pp .ge. 0.50d0) cycle + !> pi atom center sigma (nl,nm,nr) + !> from 3 atoms A-B=C, find central atom (B) + call threeoutfour(i1,i2,j1,j2,nl,nm,nr) + if (nl .eq. 0) cycle !> fall back =do nothing if assignment fails + vec1(1:3) = xyz(1:3,nm) + vec2(1:3) = (xyz(1:3,nl)+xyz(1:3,nr))*0.5 + dtot(1:3) = rklmo(1:3,pmo) + call calcrotation(dtot,vec2,vec1-vec2,pi) !> project upper pi center to other bond + k = k+1 + rklmo(1:3,n+k) = dtot(1:3) + rklmo(4:5,n+k) = rklmo(4:5,pmo) + + !> add to screen file, protomer search + if (pr_local) then +! write (iscreen,*) nat+1 +! write (iscreen,*) +! do ii = 1,nat +! write (iscreen,'(a2,3F24.10)') i2e(at(ii)),xyz(1:3,ii)*autoaa +! end do +! write (iscreen,'(a2,3F24.10)') i2e(1),dtot(1:3)*autoaa + end if + + imo = lneigh(1,pmo) + vec1(1:3) = xyz(1:3,nm) + vec2(1:3) = (xyz(1:3,nl)+xyz(1:3,nr))*0.5_wp + dtot(1:3) = rklmo(1:3,imo) + call calcrotation(dtot,vec2,vec1-vec2,pi) !> project bottom pi center to other bond + k = k+1 + rklmo(1:3,n+k) = dtot(1:3) + rklmo(4:5,n+k) = rklmo(4:5,imo) + rklmo(5,smo) = 0 ! remove sigma + + !> add to screen file, protomer search + if (pr_local) then +! write (iscreen,*) nat+1 +! write (iscreen,*) +! do ii = 1,nat +! write (iscreen,'(a2,3F24.10)') i2e(at(ii)),xyz(1:3,ii)*autoaa +! end do +! write (iscreen,'(a2,3F24.10)') i2e(1),dtot(1:3)*autoaa + end if + + m = m+1 + end do + end do + new = k + + if (pr_local) then +! close (iscreen) + end if + deallocate (wbo) + + end if + !ccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + !> collect possible protonation sites for output + nprot = 0 + k = new+n + do i=1,k + if(rklmo(5,i) > 1) nprot=nprot+1 + enddo + if(nprot > 0)then + allocate(protxyz(3,nprot), source=0.0_wp) + m=0 + do i=1,k + if(rklmo(5,i) > 1)then + m=m+1 + protxyz(1:3,m) = rklmo(1:3,i) + endif + enddo + endif + + !ccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + !> If the normal xtb mode is used with --lmo, pr_local is true and + ! a file with all information is written. + ! If the docking mode is used, information are stored intrenally + +! if (pr_local) then +! call open_file(ilmoi,set%lmoinfo_fname,'w') +! write (ilmoi,*) nat +! if (pr_local) call open_file(icent,'lmocent.coord','w') +! if (pr_local) write (icent,'(''$coord'')') +! do i = 1,nat +! write (ilmoi,'(i2,3F20.10,E18.10)') at(i),xyz(1:3,i),q(i) +! if (pr_local) write (icent,'(3F24.10,5x,a2)') xyz(1:3,i),i2e(at(i)) +! end do +! k = 0 +! do i = 1,n+new +! if (int(rklmo(5,i)) .gt. 0) k = k+1 +! end do +! write (ilmoi,'(i5,5f14.8)') k,diptot,enlumo,enhomo +! do i = 1,n+new +! if (int(rklmo(5,i)) .gt. 0) then +! write (ilmoi,'(i2,3F20.10,f14.8)') int(rklmo(5,i)),rklmo(1:3,i) +! if (pr_local) write (icent,'(3F24.10,5x,a2)') rklmo(1:3,i),i2e(2) +! end if +! end do +! write (ilmoi,'(10(F10.6))') (qhl(i,1),i=1,nat) ! HOMO atom pop +! write (ilmoi,'(10(F10.6))') (qhl(i,2),i=1,nat) ! LUMO atom pop +! if (pr_local) write (icent,'(''$end'')') +! call close_file(ilmoi) +! if (pr_local) call close_file(icent) +! end if + + if (pr_local) then +! write (*,*) 'files:' +! write (*,*) 'coordprot.0/xtbscreen.xyz/xtblmoinfo/lmocent.coord' +! write (*,*) 'with protonation site input, xtbdock and' +! write (*,*) 'LMO center info written' +! write (*,*) + end if + + deallocate (xcen,cca,d,f,qmo,ecent,rklmo) + + end subroutine local + +!========================================================================================! + subroutine lmotype(n,at,xyz,ex,ey,ez,ia1,ia2,xcen,modi,pithr,typ) +!****************************************************** +!* lmotype +!* determine type of LMO based on empirical thresholds +!****************************************************** + implicit none + integer,intent(in) :: n,ia1,ia2,at(n) + integer,intent(out) :: typ + logical,intent(in) :: modi + real(wp),intent(in) :: xyz(3,n) + real(wp),intent(inout) :: ex,ey,ez + real(wp),intent(in) :: xcen,pithr + real(wp) :: r1,r2,r,r0,f,norm + + if (xcen .lt. 1.3333333_wp) then + typ = 2 ! LP + f = -2.2_wp + if (modi) then + norm = sqrt((xyz(1,ia1)-ex)**2+(xyz(2,ia1)-ey)**2 & + & +(xyz(3,ia1)-ez)**2) + ! the LP is on the atom e.g. in TM compounds so just shift it away + if (norm .lt. 0.2) then + call shiftlp(n,at,ia1,xyz,ex,ey,ez) + else + ! put it on the line along atom...LP center + ex = ex+f*(xyz(1,ia1)-ex) + ey = ey+f*(xyz(2,ia1)-ey) + ez = ez+f*(xyz(3,ia1)-ez) + end if + end if + else + r1 = sqrt((xyz(1,ia1)-ex)**2+(xyz(2,ia1)-ey)**2+(xyz(3,ia1)-ez)**2) + r2 = sqrt((xyz(1,ia2)-ex)**2+(xyz(2,ia2)-ey)**2+(xyz(3,ia2)-ez)**2) + r = r1+r2 + r0 = sqrt(sum((xyz(:,ia1)-xyz(:,ia2))**2)) + if (r/r0 .gt. 1.04_wp) then + typ = 3 ! pi + if (xcen .gt. pithr) typ = 4 + else + typ = 1 ! sigma + end if + end if + + end subroutine lmotype + +!========================================================================================! + subroutine mocent(n,ndim,ihomo,x,s,qmo,xcen,aoat2) + implicit none + integer :: n,ndim,ihomo + integer,intent(in) :: aoat2(ndim) + real(wp) :: xcen(ihomo),x(ndim,ndim),s(ndim,ndim) + real(wp) :: qmo(n,ihomo) + + integer :: ia,ib,m,i,j,k + real(wp) :: dum,tmp + + qmo = 0 + !$OMP PARALLEL PRIVATE (m,j,ia,k,ib,tmp,dum) & + !$OMP& SHARED(qmo,xcen,X,S) & + !$OMP& DEFAULT(SHARED) + !$OMP DO + do m = 1,ihomo + do j = 1,ndim + ia = aoat2(j) + do k = 1,j-1 + ib = aoat2(k) + tmp = X(j,m)*X(k,m)*s(k,j) + qmo(ia,m) = qmo(ia,m)+tmp + qmo(ib,m) = qmo(ib,m)+tmp + end do + qmo(ia,m) = qmo(ia,m)+X(j,m)*X(j,m)*s(j,j) + end do + dum = 0 + do i = 1,n + dum = dum+qmo(i,m)**2 + end do + xcen(m) = 1.0/(1.0d-8+dum) + end do + !$OMP END DO + !$OMP END PARALLEL + + end subroutine mocent + + SUBROUTINE lmosort(ncent,ihomo,imo,imem,qmo) + integer :: ncent,ihomo,imo + integer :: imem(ncent) + real(wp) :: qmo(ncent,ihomo) + integer :: ii,i,k,j,ihilf + real(wp) :: pp + do ii = 2,ncent + i = ii-1 + k = i + pp = qmo(i,imo) + do j = ii,ncent + if (qmo(j,imo) .lt. pp) cycle + k = j + pp = qmo(j,imo) + end do + if (k .eq. i) cycle + qmo(k,imo) = qmo(i,imo) + qmo(i,imo) = pp + + ihilf = imem(i) + imem(i) = imem(k) + imem(k) = ihilf + end do + + end subroutine lmosort + +!========================================================================================! + SUBROUTINE lmosort2(n,eps,d,ecent) + real(wp) :: d(n,n),eps(n),ecent(n,3) + integer :: n + integer :: ii,j,i,k + real(wp) :: hilf,pp + do ii = 2,n + i = ii-1 + k = i + pp = eps(i) + do j = ii,n + if (eps(j) .gt. pp) cycle + k = j + pp = eps(j) + end do + if (k .eq. i) cycle + eps(k) = eps(i) + eps(i) = pp + + do j = 1,n + hilf = d(j,i) + d(j,i) = d(j,k) + d(j,k) = hilf + end do + + do j = 1,3 + hilf = ecent(i,j) + ecent(i,j) = ecent(k,j) + ecent(k,j) = hilf + end do + + end do + + end subroutine lmosort2 + +!========================================================================================! + subroutine lmoneigh(n,rk,ecent,aneigh,neigh) + implicit none + integer :: n,neigh(4,n),aneigh(2,*) + real(wp) :: ecent(n,3) + real(sp) :: rk(5,2*n) + + real(wp) rr(n) + integer :: i,j,ind(n),i1,i2,j1,j2 + external qsort + + do i = 1,n + i1 = aneigh(1,i) + i2 = aneigh(2,i) + do j = 1,n + j1 = aneigh(1,j) + j2 = aneigh(2,j) + rr(j) = (ecent(i,1)-ecent(j,1))**2 & + & +(ecent(i,2)-ecent(j,2))**2 & + & +(ecent(i,3)-ecent(j,3))**2 + if (int(rk(5,j)) .eq. 2) rr(j) = 1.d+42 ! LP + if (i1 .eq. j1.and.i2 .eq. j2) rr(j) = -1.0 ! same bond? + if (i1 .eq. j2.and.i2 .eq. j1) rr(j) = -1.0 + if (i .eq. j) rr(j) = 2.d+42 + ind(j) = j + end do + call qsort(rr,1,n,ind) + neigh(1,i) = ind(1) + neigh(2,i) = ind(2) + neigh(3,i) = ind(3) + neigh(4,i) = ind(4) + end do + + end subroutine lmoneigh + +!========================================================================================! + subroutine atomneigh(n,nat,xyz,ecent,neigh) + implicit none + integer,intent(in) :: n,nat + integer,intent(inout) :: neigh(2,n) + real(wp),intent(in) :: ecent(n,3) + real(wp),intent(in) :: xyz(3,nat) + + real(wp) :: rr(nat) + integer :: i,j,ind(nat) + external qsort + do i = 1,n + do j = 1,nat + rr(j) = (ecent(i,1)-xyz(1,j))**2 & + & +(ecent(i,2)-xyz(2,j))**2 & + & +(ecent(i,3)-xyz(3,j))**2 + ind(j) = j + end do + call qsort(rr,1,nat,ind) + neigh(1,i) = ind(1) + neigh(2,i) = ind(2) + end do + + end subroutine atomneigh + +!========================================================================================! + pure function bndcheck(nat,list,i1,i2) result(check) + integer,intent(in) :: nat,list(nat),i1,i2 + logical :: check + check = .false. + if (list(i1) .eq. 1.and.list(i2) .eq. 1) check = .true. + end function bndcheck + +!========================================================================================! + subroutine irand3(n1,n2,n3) + integer,intent(out) :: n1,n2,n3 + integer :: irand + real(sp) :: x + + call random_number(x) + irand = int(3.1*x) + if (irand .lt. 1) irand = 1 + if (irand .gt. 3) irand = 3 + n1 = irand +10 call random_number(x) + irand = int(3.1*x) + if (irand .lt. 1) irand = 1 + if (irand .gt. 3) irand = 3 + if (irand .ne. n1) then + n2 = irand + else + goto 10 + end if +20 call random_number(x) + irand = int(3.1*x) + if (irand .lt. 1) irand = 1 + if (irand .gt. 3) irand = 3 + if (irand .ne. n1.and.irand .ne. n2) then + n3 = irand + else + goto 20 + end if + + end subroutine irand3 + +!========================================================================================! + pure subroutine threeoutfour(i1,i2,j1,j2,n1,n2,n3) + implicit none + integer,intent(in) :: i1,i2,j1,j2 + integer,intent(out) :: n1,n2,n3 + + n1 = 0 + if (i1 .eq. j1) then + n1 = i2 + n2 = i1 + n3 = j2 + return + end if + if (i1 .eq. j2) then + n1 = i2 + n2 = i1 + n3 = j1 + return + end if + if (i2 .eq. j1) then + n1 = i1 + n2 = i2 + n3 = j2 + return + end if + if (i2 .eq. j2) then + n1 = i1 + n2 = i2 + n3 = j1 + return + end if + + end subroutine threeoutfour + +!========================================================================================! + pure subroutine calcrotation(x,ori,vec,phi) + implicit none + integer :: i,j + real(wp),intent(inout) :: x(3) + real(wp),intent(in) :: ori(3) + real(wp),intent(in) :: vec(3) + real(wp),intent(in) :: phi + real(wp) :: d(3) + real(wp) :: xtmp(3) + real(wp) :: absd + + d(1:3) = vec(1:3) + + xtmp(1:3) = x(1:3)-ori(1:3) + + absd = sqrt(d(1)**2+d(2)**2+d(3)**2) + + d(1:3) = d(1:3)/absd + + x(1) = ((d(2)**2+d(3)**2)*cos(phi)+d(1)**2)*xtmp(1)+ & + & (d(1)*d(2)*(1-cos(phi))-d(3)*sin(phi))*xtmp(2)+ & + & (d(1)*d(3)*(1-cos(phi))+d(2)*sin(phi))*xtmp(3) + x(2) = (d(1)*d(2)*(1-cos(phi))+d(3)*sin(phi))*xtmp(1)+ & + & ((d(1)**2+d(3)**2)*cos(phi)+d(2)**2)*xtmp(2)+ & + & (d(2)*d(3)*(1-cos(phi))-d(1)*sin(phi))*xtmp(3) + x(3) = (d(1)*d(3)*(1-cos(phi))-d(2)*sin(phi))*xtmp(1)+ & + & (d(2)*d(3)*(1-cos(phi))+d(1)*sin(phi))*xtmp(2)+ & + & ((d(1)**2+d(2)**2)*cos(phi)+d(3)**2)*xtmp(3) + x(1:3) = x(1:3)+ori(1:3) + + end subroutine calcrotation + +!========================================================================================! + pure subroutine piorient(a,b,flip) + implicit none + real(wp),intent(in) :: a(3),b(3) + logical,intent(out) :: flip + integer :: i,j + real(wp) :: d(3),d2(3) + flip = .false. + d = 0.0d0 + j = 0 + + do i = 1,3 + d(i) = a(i)-b(i) + d2(i) = d(i)*d(i) + end do + j = maxloc(d2,1) + if (d(j) .lt. 0.0d0) flip = .true. + return + end subroutine piorient + +!========================================================================================! + + subroutine Dints(n,nbf,xyz,S1,S2,S3, & + & aoat,nprim,alp,lao,cont) + use xtb_intpack,only:opab1,propa + implicit none + !> INPUT + integer,intent(in) :: n + integer,intent(in) :: nbf + real(wp),intent(in) :: xyz(3,n) + integer,intent(in) :: aoat(nbf) !> mapping BF -> atom + integer,intent(in) :: nprim(nbf) !> mappinf BF -> primitive number + real(wp),intent(in) :: alp(9*nbf) !> mapping primitive -> primitive exponent + integer,intent(in) :: lao(nbf) !> mapping BF -> azimudal quantum number l + real(wp),intent(in) :: cont(9*nbf) !> mapping primitive -> contraction coeffient + !> OUTPUT + real(wp),intent(out) :: S1(nbf*(nbf+1)/2) + real(wp),intent(out) :: S2(nbf*(nbf+1)/2) + real(wp),intent(out) :: S3(nbf*(nbf+1)/2) + !> LOCAL + integer :: i,j,k,l + integer :: iprimcount,jprimcount + integer :: npri,nprj + integer :: ii,iii,jj,jjj,ll,m,li,lj,mm,nn + real(wp) :: xyza(3),xyzb(3),rab,est,ss,sss,gama,arg + real(wp) :: point(3),gm2,ttt(3),tt1,tt2,tt3,intcut + + intcut = 20.0_wp + + point = 0.0_wp + s1 = 0.0_wp + s2 = 0.0_wp + s3 = 0.0_wp + + k = 0 + iprimcount = 0 + do i = 1,nbf + !> aufpunkt i + xyza(1:3) = xyz(1:3,aoat(i)) + !> #prims + npri = nprim(i) + jprimcount = 0 + do j = 1,i + k = k+1 + nprj = nprim(j) + !> aufpunkt j + xyzb(1:3) = xyz(1:3,aoat(j)) + rab = sum((xyza-xyzb)**2) + if (rab .gt. 200) goto 42 !> cut-off gives crambin dipole accurate to 1d-3 Deb + !> prim loop + tt1 = 0.0_wp + tt2 = 0.0_wp + tt3 = 0.0_wp + do ii = 1,npri + iii = iprimcount+ii + do jj = 1,nprj + jjj = jprimcount+jj + gama = 1.0_wp/(alp(iii)+alp(jjj)) + est = rab*alp(iii)*alp(jjj)*gama + !> cutoff + if (est .lt. intcut) then + ttt = 0 + call propa(opab1,xyza,xyzb,point,alp(iii),alp(jjj),& + & lao(i),lao(j),ttt,3) + tt1 = tt1+ttt(1)*cont(iii)*cont(jjj) + tt2 = tt2+ttt(2)*cont(iii)*cont(jjj) + tt3 = tt3+ttt(3)*cont(iii)*cont(jjj) + end if + end do + end do + s1(k) = tt1 + s2(k) = tt2 + s3(k) = tt3 +42 jprimcount = jprimcount+nprj + end do + iprimcount = iprimcount+npri + end do + + end subroutine Dints + +!========================================================================================! + + subroutine onetri(ity,s,s1,array,n,ival) + ! ******designed for abelian groups only****** + ! + ! calling sequence: + ! ity =1 sym ao + ! =-1 anti sym ao + ! s input property matrix over ao's + ! s1 transformed integrals on output + ! s2 scratch arrays large enough to hold a square matrix + ! array mo matrix over so's + ! n linear dimension of arrays + ! bernd hess, university of bonn, january 1991 + integer,intent(in) :: n + integer,intent(in) :: ity + integer,intent(in) :: ival + real(wp),intent(in) :: s(*) + real(wp),intent(inout) :: s1(*) + real(wp),intent(in) :: array(n,ival) + real(wp) :: s2(n,n) + external :: dsymm,dgemm ! can't use wrappers due to change of leading dim + + ! + ! determine if we have an antisymmetric integral + if (ity /= -1) then + ! + ! blow up symmetric matrix s + call blowsy(ity,s,s1,n) + + ! + ! transformation of s + call dsymm('l','l',n,ival,1.d0,s1,n,array,n,0.d0,s2,n) + call dgemm('t','n',ival,ival,n,1.d0,array,n,s2,n,0.d0,s1,ival) + else + ! + ! blow up anti-symmetric matrix s + call blowsy(ity,s,s1,n) + ! + ! transformation of s + call dgemm('n','n',n,ival,n,1.d0,s1,n,array,n,0.d0,s2,n) + call dgemm('t','n',ival,ival,n,1.d0,array,n,s2,n,0.d0,s1,ival) + end if + end subroutine + +!========================================================================================! + + pure subroutine blowsy(ity,a,b,n) + implicit none + ! blow up symmetric or antisymmetric matrix to full size + integer,intent(in) :: ity + integer,intent(in) :: n + real(wp),intent(in) :: a(*) + real(wp),intent(out) :: b(n,n) + integer :: ij,i,j + ! determine if we have an antisymmetric integral + + ij = 0 + if (ity .eq. -1) then + do i = 1,n + do j = 1,i-1 + ij = ij+1 + b(j,i) = -a(ij) + b(i,j) = a(ij) + end do + ij = ij+1 + b(i,i) = 0.d0 + end do + else + do i = 1,n + do j = 1,i-1 + ij = ij+1 + b(j,i) = a(ij) + b(i,j) = a(ij) + end do + ij = ij+1 + b(i,i) = a(ij) + end do + end if + end subroutine blowsy + +!========================================================================================! + + subroutine cao2saop(nbf,nao,s,lao) + implicit none + !> INPUT + integer,intent(in) :: nbf !> number of BF + integer,intent(inout) :: nao !> number of MOs + integer,intent(in) :: lao(nbf) !> mapping BF -> azimudal quantum number l + !> OUTPUT + real(wp),intent(inout) :: s(nbf*(nbf+1)/2) + !> LOCAL + real(wp) :: sspher + integer :: lll(20),firstd(nbf) + integer :: i,j,k,ii,jj,kk,iii,jjj,li,lj,mm,nn,m,n + data lll/1,2,2,2,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4/ + + real(wp) :: trafo(6,6) + real(wp),allocatable :: sneu(:) + + trafo = 0.0d0 + ! dS + trafo(1,1) = 1.d0/sqrt(3.d0)*sqrt(3.d0/5.d0) + trafo(2,1) = 1.d0/sqrt(3.d0)*sqrt(3.d0/5.d0) + trafo(3,1) = 1.d0/sqrt(3.d0)*sqrt(3.d0/5.d0) + ! dx2-y2 + trafo(1,2) = 1.d0/sqrt(2.d0)*sqrt(3.d0/2.d0) + trafo(2,2) = -1.d0/sqrt(2.d0)*sqrt(3.d0/2.d0) + ! dz2 + trafo(1,3) = 0.50d0 + trafo(2,3) = 0.50d0 + trafo(3,3) = -1.0d0 + ! rest + trafo(4,4) = 1.0d0 + trafo(5,5) = 1.0d0 + trafo(6,6) = 1.0d0 + + nao = 0 + firstd = 0 + i = 1 +42 if (lao(i) .gt. 4.and.lao(i) .le. 10) then + nao = nao+1 + firstd(i:i+5) = i + i = i+5 + end if + i = i+1 + if (i .lt. nbf) goto 42 + + if (nao .eq. 0) then + nao = nbf + return + end if + + allocate (sneu(nbf*(nbf+1)/2)) + + sneu = s + + k = 0 + do i = 1,nbf + li = lll(lao(i)) + do j = 1,i + lj = lll(lao(j)) + k = k+1 + ! d-d + if (li .eq. 3.and.lj .eq. 3) then + ii = lao(i)-4 + jj = lao(j)-4 + sspher = 0 + do m = 1,6 + mm = firstd(i)-1+m + do n = 1,6 + nn = firstd(j)-1+n + sspher = sspher+trafo(m,ii)*trafo(n,jj)*s(lin(mm,nn)) + end do + end do + sneu(k) = sspher + end if + ! d-sp + if (li .eq. 3.and.lj .le. 2) then + ii = lao(i)-4 + sspher = 0 + do m = 1,6 + mm = firstd(i)-1+m + sspher = sspher+trafo(m,ii)*s(lin(mm,j)) + end do + sneu(k) = sspher + end if + ! sp-d + if (li .le. 2.and.lj .eq. 3) then + jj = lao(j)-4 + sspher = 0 + do n = 1,6 + nn = firstd(j)-1+n + sspher = sspher+trafo(n,jj)*s(lin(i,nn)) + end do + sneu(k) = sspher + end if + + end do + end do + + s(1:nao*(nao+1)/2) = 0 + + k = 0 + iii = 0 + do i = 1,nbf + if (lao(i) .ne. 5) iii = iii+1 + jjj = 0 + do j = 1,i + if (lao(j) .ne. 5) jjj = jjj+1 + k = k+1 + if (lao(i) .eq. 5.or.lao(j) .eq. 5) cycle + s(lin(iii,jjj)) = sneu(k) + end do + end do + + nao = nbf-nao + + deallocate (sneu) + end subroutine cao2saop + +!========================================================================================! + + integer function lin(i1,i2) + integer :: i1,i2,idum1,idum2 + idum1 = max(i1,i2) + idum2 = min(i1,i2) + lin = idum2+idum1*(idum1-1)/2 + return + end function lin + +!========================================================================================! + + subroutine shiftlp(n,at,ia1,xyz,ex,ey,ez) +! shift the LP=protonation position to an "empty" region + implicit none + integer :: n,at(n),ia1 + real(wp) :: xyz(3,n),ex,ey,ez + + real(wp) :: r0,r1,r2,exmin,eymin,ezmin,rmin,xsave,ysave,zsave,rmax + real(wp) :: x,f + integer :: icount,j + + r0 = 2.5d0 + rmin = 1.d+42 + rmax = 0 + icount = 0 + +10 f = 1.0 + call random_number(x) + if (x .lt. 0.5) f = -1.0 + call random_number(x) + ex = xyz(1,ia1)+x*f*r0 + call random_number(x) + ey = xyz(2,ia1)+x*f*r0 + call random_number(x) + ez = xyz(3,ia1)+x*f*r0 + r1 = sqrt((xyz(1,ia1)-ex)**2+(xyz(2,ia1)-ey)**2+(xyz(3,ia1)-ez)**2) + if (abs(r0-r1) .gt. 0.1) goto 10 + rmin = 1.d+42 + do j = 1,n + r2 = sqrt((xyz(1,j)-ex)**2+(xyz(2,j)-ey)**2+(xyz(3,j)-ez)**2) + if (r2 .lt. rmin.and.j .ne. ia1) then + rmin = r2 + exmin = ex + eymin = ey + ezmin = ez + end if + end do + if (rmin .gt. rmax) then + rmax = rmin + xsave = exmin + ysave = eymin + zsave = ezmin + end if + icount = icount+1 + if (icount .lt. 100) goto 10 + + ex = xsave + ey = ysave + ez = zsave + + end subroutine shiftlp + +!========================================================================================! +!========================================================================================! +end module mo_localize diff --git a/src/qmhelpers/lopt.f90 b/src/qmhelpers/lopt.f90 new file mode 100644 index 00000000..5119d9c3 --- /dev/null +++ b/src/qmhelpers/lopt.f90 @@ -0,0 +1,233 @@ +! Copyright (C) 2017-2020 Stefan Grimme, Sebastian Ehlert (xtb) +! Copyright (C) 2024 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 +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! crest is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with crest. If not, see . +! +! Routines were adapted from the xtb code (github.com/grimme-lab/xtb) +! under the Open-source software LGPL-3.0 Licencse. +!================================================================================! + +module lopt_mod + use iso_fortran_env,only:wp => real64 + public + + logical,private,parameter :: debug = .false. + +!================================================================================! +!================================================================================! +contains !> MODULE PROCEDURES START HERE +!================================================================================! +!================================================================================! + + subroutine lopt(init,n,no,accr,op,d) +!******************************************************************** +!* Original subroutine for Boys localization by R.Ahlrichs 2/85 +!* +!* to determine vectors |i>, i=1,n +!* which maximize +!* f = Σ(i,j=1,n)Σ(k=1,no) (-)**2 +!* equivalent to maximization of +!* g = Σ(i=1,n)Σ(k=1,no) **2 + **2 +!* equivalent to minimization of +!* h = Σ(i**2 +!* stationarity condition +!* Σ(k) op(ij,k)*(op(ii,k)-op(jj,k)) = 0 all i.ne.j +!* +!* input +!* n = dimension of lin. space +!* no = # of operators in f or g +!* op(ij,k) matrix representation of k'th operator +!* symmetric and stored triangular,nn1=n(n+1)/2 +!* thi = rotation threshold (1.-6 to 1.-12 sngl or dbl prc) +!* nboys = max. nb. of sweeps +!* result +!* d - contains the solution vectors as columns +!* op - the transformed operators +!* +!******************************************************************* + implicit real(wp) (a-h,o-z) + !> INPUT + logical,intent(in) :: init + integer,intent(in) :: n,no + real(wp),intent(in) :: accr + !> OUTPUT + real(wp),intent(inout) :: op(n*(n+1)/2,no) + real(wp),intent(out) :: d(n,n) + !> LOCAL + integer :: auxi + real(wp) :: auxvec(n) + + a0 = 0.d0 + a1 = 1.d0 + a2 = 2.d0 + a10 = 10.d0 + aqu = 0.25d0 + d(:,:) = 0.0_wp +!>--- accr=1.d-6 ! originally 10-14 but this setting saves a factor of 3 ! + thi = 1.d-10 + th = dmax1(thi,accr) + + first = 0 + +!>--- set d to unit matrix + if (init) then + if (debug) write (*,*) 'initialization of trafo matrix to unity' + do i = 1,n + do j = 1,n +10 d(j,i) = a0-accr*i+accr*i*j + end do +20 d(i,i) = a1 + end do + end if + + if (n .eq. 1) return + + nboys = 20000 + nn1 = n*(n+1)/2 + +!>--- opn = operator norm needed for cut off threshold + opn = a0 + do k = 1,no + do i = 1,nn1 +21 opn = opn+op(i,k)**2 + end do +22 continue + end do + if (opn .eq. a0) return + opn = dsqrt(opn/dble(nn1*no))*accr + + ths = aqu + it = 0 + do isw = 1,nboys +!>--- loop over sweeps + thsw = a0 + ii = 1 + ij = 2 + do i = 2,n + ii = ii+i + jj = 0 + im1 = i-1 + ip1 = i+1 + li = ii-i + do j = 1,im1 +!>--- loop over i,j rotations +!> maximize sum(k) **2 + **2 +!> like in jacobi + jj = jj+j +!>--- get ax , b determining the tan of the rotation angle + ax = a0 + b = a0 + do k = 1,no + aux1 = op(jj,k)-op(ii,k) + aux2 = op(ij,k) + ax = ax+aux1*aux2 +30 b = b+aqu*aux1**2-aux2**2 + end do +!>--- rotation angle x0 + if (dabs(ax)+dabs(b) .lt. opn) go to 100 + x0 = aqu*datan2(ax,b) + thsw = dmax1(thsw,dabs(x0)) + if (dabs(x0) .lt. ths) go to 100 +!>--- rotation + it = it+1 + s = dsin(x0) + s2 = s*s + c2 = a1-s2 + c = dsqrt(c2) + do l = 1,n + aux1 = d(l,j) + d(l,j) = c*d(l,j)+s*d(l,i) + d(l,i) = c*d(l,i)-s*aux1 +40 continue + end do + jm1 = j-1 + cs = c*s + c2ms2 = c2-s2 + lj = jj-j + jp1 = j+1 + + ! !$OMP PARALLEL PRIVATE (k,l,aux1,aux2,il,jl) & + ! !$OMP& SHARED(op) & + ! !$OMP& DEFAULT(SHARED) + ! !$OMP DO + do k = 1,no + aux1 = op(ii,k) + aux2 = op(jj,k) + op(ii,k) = aux1*c2+aux2*s2-a2*op(ij,k)*cs + op(jj,k) = aux1*s2+aux2*c2+a2*op(ij,k)*cs + op(ij,k) = cs*(aux1-aux2)+op(ij,k)*c2ms2 + if (jm1 .le. 0) go to 60 + + ! do l=1,jm1 + ! aux1=c*op(li+l,k)-s*op(lj+l,k) + ! op(lj+l,k)=s*op(li+l,k)+c*op(lj+l,k) + ! op(li+l,k)=aux1 + ! end do + auxvec(1:jm1) = c*op(li+1:li+jm1,k)-s*op(lj+1:lj+jm1,k) + op(lj+1:lj+jm1,k) = op(lj+1:lj+jm1,k)*c + op(lj+1:lj+jm1,k) = op(lj+1:lj+jm1,k)+s*op(li+1:li+jm1,k) + op(li+1:li+jm1,k) = auxvec(1:jm1) + +60 if (im1 .eq. j) go to 80 + jl = jj+j + + do l = jp1,im1 + aux1 = c*op(li+l,k)-s*op(jl,k) + op(jl,k) = c*op(jl,k)+s*op(li+l,k) + op(li+l,k) = aux1 +70 jl = jl+l + end do + +80 if (i .eq. n) go to 95 + jl = ii+j + il = ii+i + + do l = ip1,n + aux1 = c*op(jl,k)+s*op(il,k) + op(il,k) = c*op(il,k)-s*op(jl,k) + op(jl,k) = aux1 + il = il+l +90 jl = jl+l + end do + +95 continue + end do + ! !$OMP END DO + ! !$OMP END PARALLEL + +100 ij = ij+1 + end do +110 ij = ij+1 + end do + +!>--- check convergence + if (thsw .le. th) go to 210 + ths = dmin1(thsw**2,ths*aqu) + ths = dmax1(th,ths) + if (mod(isw,50) .eq. 0.AND.debug)& + & write (*,'(''iteration'',i7,2x,''convergence'',d16.8)') isw,thsw +200 continue + end do + if (debug) write (*,300) isw,thsw + return +210 if (debug) write (*,310) isw,thsw + + return +300 format(/' not converged in',i7,' iterations, threshold :',d16.8) +310 format(/' converged in',i7,' iterations, threshold : ',d16.8) + end subroutine lopt + +!========================================================================================! +!========================================================================================! +end module lopt_mod diff --git a/src/qmhelpers/meson.build b/src/qmhelpers/meson.build index 0bb843b0..0716e307 100644 --- a/src/qmhelpers/meson.build +++ b/src/qmhelpers/meson.build @@ -16,4 +16,7 @@ srcs += files( 'wiberg_mayer.f90', + 'intpack.f90', + 'lopt.f90', + 'local.f90', ) diff --git a/src/sortens.f90 b/src/sortens.f90 index 83cd384e..bcc2cc3c 100644 --- a/src/sortens.f90 +++ b/src/sortens.f90 @@ -52,7 +52,7 @@ subroutine sort_ens(sort,infile,verbose) ochan = 6 !stdout - !sort=env%ptb + !sort=env%protb associate (popthr => sort%popthr, & & allowFrag => sort%allowFrag,ewin => sort%ewin, & diff --git a/src/strucreader.f90 b/src/strucreader.f90 index 07cfb5a1..56078d52 100644 --- a/src/strucreader.f90 +++ b/src/strucreader.f90 @@ -205,6 +205,7 @@ module strucrd procedure :: dihedral => coord_getdihedral !> calculate dihedral angle between four atoms procedure :: cutout => coord_getcutout !> create a substructure procedure :: get_CN => coord_get_CN !> calculate coordination number + procedure :: get_z => coord_get_z !> calculate nuclear charge end type coord !=========================================================================================! !ensemble class. contains all structures of an ensemble @@ -1446,6 +1447,22 @@ subroutine coord_get_CN(self,cn,cn_type,cn_thr,dcndr) & cntype=cn_type, cnthr=cn_thr, dcndr=dcndr) end subroutine coord_get_CN +!==================================================================! + + subroutine coord_get_z(self,z) + implicit none + class(coord) :: self + real(wp),intent(out),allocatable :: z(:) + integer :: i,j,k + if(self%nat <= 0) return + if(.not.allocated(self%xyz) .or. .not.allocated(self%at)) return + allocate(z(self%nat), source=0.0_wp) + do i=1,self%nat + z(i) = real(self%at(i),wp) - real(ncore(self%at(i))) + if (self%at(i) > 57 .and. self%at(i) < 72) z(i) = 3.0_wp + enddo + end subroutine coord_get_z + !=========================================================================================! !=========================================================================================! ! 3. ROUTINES FOR WRITING STRUCTURES AND CONVERTING THEM @@ -2011,6 +2028,32 @@ function convertlable(s) call move_alloc(sout,convertlable) end function convertlable +!=============================================================! +pure elemental integer function ncore(at) + integer,intent(in) :: at + if(at.le.2)then + ncore=0 + elseif(at.le.10)then + ncore=2 + elseif(at.le.18)then + ncore=10 + elseif(at.le.29)then !zn + ncore=18 + elseif(at.le.36)then + ncore=28 + elseif(at.le.47)then + ncore=36 + elseif(at.le.54)then + ncore=46 + elseif(at.le.71)then + ncore=54 + elseif(at.le.79)then + ncore=68 + elseif(at.le.86)then + ncore=78 + endif +end function ncore + !============================================================! ! e2i is used to map the element (as a string) to integer !============================================================! From 44a9ff61cdd0b7022f1ab518569ae01cc53d83e9 Mon Sep 17 00:00:00 2001 From: Philipp Pracht Date: Fri, 19 Jul 2024 02:01:32 +0200 Subject: [PATCH 2/5] added interfaces for parallel optimization driver --- src/algos/optimization.f90 | 1 + src/algos/parallel.f90 | 195 +++++++++++++++++++++----------- src/algos/protonate.f90 | 1 + src/algos/refine.f90 | 1 + src/algos/search_1.f90 | 1 + src/algos/search_conformers.f90 | 1 + src/algos/search_mecp.f90 | 1 + src/qmhelpers/local.f90 | 13 ++- 8 files changed, 141 insertions(+), 73 deletions(-) diff --git a/src/algos/optimization.f90 b/src/algos/optimization.f90 index 793e1595..cb885348 100644 --- a/src/algos/optimization.f90 +++ b/src/algos/optimization.f90 @@ -140,6 +140,7 @@ subroutine crest_ensemble_optimization(env,tim) use crest_calculator use strucrd use optimize_module + use parallel_interface implicit none type(systemdata),intent(inout) :: env type(timer),intent(inout) :: tim diff --git a/src/algos/parallel.f90 b/src/algos/parallel.f90 index 45b1e10f..13fd6f10 100644 --- a/src/algos/parallel.f90 +++ b/src/algos/parallel.f90 @@ -19,6 +19,57 @@ !> A collection of routines to set up OMP-parallel runs of MDs and optimizations. +!========================================================================================! +!========================================================================================! +!> Interfaces to handle optional arguments +!========================================================================================! +!========================================================================================! +module parallel_interface +!******************************************************* +!* module to load an interface to the parallel routines +!* mandatory to handle any optional input arguments +!******************************************************* + implicit none + interface + subroutine crest_sploop(env,nat,nall,at,xyz,eread) + use crest_parameters,only:wp,stdout,sep + use crest_calculator + use omp_lib + use crest_data + use strucrd + use optimize_module + use iomod,only:makedir,directory_exist,remove + implicit none + type(systemdata),intent(inout) :: env + real(wp),intent(inout) :: xyz(3,nat,nall) + integer,intent(in) :: at(nat) + real(wp),intent(inout) :: eread(nall) + integer,intent(in) :: nat,nall + end subroutine crest_sploop + end interface + + interface + subroutine crest_oloop(env,nat,nall,at,xyz,eread,dump,customcalc) + use crest_parameters,only:wp,stdout,sep + use crest_calculator + use omp_lib + use crest_data + use strucrd + use optimize_module + use iomod,only:makedir,directory_exist,remove + use crest_restartlog,only:trackrestart,restart_write_dummy + implicit none + type(systemdata),target,intent(inout) :: env + real(wp),intent(inout) :: xyz(3,nat,nall) + integer,intent(in) :: at(nat) + real(wp),intent(inout) :: eread(nall) + integer,intent(in) :: nat,nall + logical,intent(in) :: dump + type(calcdata),intent(in),target,optional :: customcalc + end subroutine crest_oloop + end interface +end module parallel_interface + !========================================================================================! !========================================================================================! !> Routines for concurrent singlepoint evaluations @@ -93,7 +144,7 @@ subroutine crest_sploop(env,nat,nall,at,xyz,eread) call crest_oloop_pr_progress(env,nall,0) !>--- shared variables - allocate (grads(3,nat,T), source=0.0_wp) + allocate (grads(3,nat,T),source=0.0_wp) c = 0 !> counter of successfull optimizations k = 0 !> counter of total optimization (fail+success) z = 0 !> counter to perform optimization in right order (1...nall) @@ -151,15 +202,15 @@ subroutine crest_sploop(env,nat,nall,at,xyz,eread) !>--- prepare some summary printout percent = float(c)/float(nall)*100.0_wp - write(atmp,'(f5.1,a)') percent,'% success)' + write (atmp,'(f5.1,a)') percent,'% success)' write (stdout,'(">",1x,i0,a,i0,a,a)') c,' of ',nall,' structures successfully evaluated (', & & trim(adjustl(atmp)) write (atmp,'(">",1x,a,i0,a)') 'Total runtime for ',nall,' singlepoint calculations:' call profiler%write_timing(stdout,1,trim(atmp),.true.) runtime = profiler%get(1) - write(atmp,'(f16.3,a)') runtime/real(nall,wp),' sec' - write(stdout,'(a,a,a)') '> Corresponding to approximately ',trim(adjustl(atmp)), & - & ' per processed structure' + write (atmp,'(f16.3,a)') runtime/real(nall,wp),' sec' + write (stdout,'(a,a,a)') '> Corresponding to approximately ',trim(adjustl(atmp)), & + & ' per processed structure' deallocate (grads) call profiler%clear() @@ -173,29 +224,33 @@ end subroutine crest_sploop !> Routines for concurrent geometry optimization !========================================================================================! !========================================================================================! -subroutine crest_oloop(env,nat,nall,at,xyz,eread,dump) +subroutine crest_oloop(env,nat,nall,at,xyz,eread,dump,customcalc) !*************************************************************** !* subroutine crest_oloop !* This subroutine performs concurrent geometry optimizations !* for the given ensemble. Inputs xyz and eread are overwritten +!* env - contains parallelization and other program settings +!* dump - decides on whether to dump an ensemble file +!* customcalc - customized (optional) calculation level data !* !* IMPORTANT: xyz should be in Bohr(!) for this routine !*************************************************************** - use crest_parameters,only:wp,stdout,sep + use crest_parameters,only:wp,stdout,sep use crest_calculator use omp_lib use crest_data use strucrd use optimize_module use iomod,only:makedir,directory_exist,remove - use crest_restartlog, only: trackrestart,restart_write_dummy + use crest_restartlog,only:trackrestart,restart_write_dummy implicit none - type(systemdata),intent(inout) :: env + type(systemdata),target,intent(inout) :: env real(wp),intent(inout) :: xyz(3,nat,nall) integer,intent(in) :: at(nat) real(wp),intent(inout) :: eread(nall) integer,intent(in) :: nat,nall logical,intent(in) :: dump + type(calcdata),intent(in),target,optional :: customcalc type(coord),allocatable :: mols(:) type(coord),allocatable :: molsnew(:) @@ -208,35 +263,42 @@ subroutine crest_oloop(env,nat,nall,at,xyz,eread,dump) integer :: thread_id,vz,job character(len=80) :: atmp real(wp) :: percent,runtime - + type(calcdata),pointer :: mycalc type(timer) :: profiler !>--- decide wether to skip this call - if(trackrestart(env))then - call restart_write_dummy(ensemblefile) - return - endif + if (trackrestart(env)) then + call restart_write_dummy(ensemblefile) + return + end if + +!>--- check which calc to use + if(present(customcalc))then + mycalc => customcalc + else + mycalc => env%calc + endif !>--- check if we have any calculation settings allocated - if (env%calc%ncalculations < 1) then + if (mycalc%ncalculations < 1) then write (stdout,*) 'no calculations allocated' return end if !>--- prepare objects for parallelization T = env%threads - allocate (calculations(T),source=env%calc) + allocate (calculations(T),source=mycalc) allocate (mols(T),molsnew(T)) do i = 1,T - do j = 1,env%calc%ncalculations - calculations(i)%calcs(j) = env%calc%calcs(j) + do j = 1,mycalc%ncalculations + calculations(i)%calcs(j) = mycalc%calcs(j) !>--- directories - ex = directory_exist(env%calc%calcs(j)%calcspace) + ex = directory_exist(mycalc%calcs(j)%calcspace) if (.not.ex) then - io = makedir(trim(env%calc%calcs(j)%calcspace)) + io = makedir(trim(mycalc%calcs(j)%calcspace)) end if write (atmp,'(a,"_",i0)') sep,i - calculations(i)%calcs(j)%calcspace = env%calc%calcs(j)%calcspace//trim(atmp) + calculations(i)%calcs(j)%calcspace = mycalc%calcs(j)%calcspace//trim(atmp) call calculations(i)%calcs(j)%printid(i,j) end do calculations(i)%pr_energies = .false. @@ -293,7 +355,7 @@ subroutine crest_oloop(env,nat,nall,at,xyz,eread,dump) !>-- geometry optimization call optimize_geometry(mols(job),molsnew(job),calculations(job),energy,grads(:,:,job),pr,wr,io) - + !$omp critical if (io == 0) then !>--- successful optimization (io==0) @@ -328,15 +390,15 @@ subroutine crest_oloop(env,nat,nall,at,xyz,eread,dump) !>--- prepare some summary printout percent = float(c)/float(nall)*100.0_wp - write(atmp,'(f5.1,a)') percent,'% success)' + write (atmp,'(f5.1,a)') percent,'% success)' write (stdout,'(">",1x,i0,a,i0,a,a)') c,' of ',nall,' structures successfully optimized (', & & trim(adjustl(atmp)) write (atmp,'(">",1x,a,i0,a)') 'Total runtime for ',nall,' optimizations:' call profiler%write_timing(stdout,1,trim(atmp),.true.) runtime = profiler%get(1) - write(atmp,'(f16.3,a)') runtime/real(nall,wp),' sec' - write(stdout,'(a,a,a)') '> Corresponding to approximately ',trim(adjustl(atmp)), & - & ' per processed structure' + write (atmp,'(f16.3,a)') runtime/real(nall,wp),' sec' + write (stdout,'(a,a,a)') '> Corresponding to approximately ',trim(adjustl(atmp)), & + & ' per processed structure' !>--- close files (if they are open) if (dump) then @@ -359,9 +421,9 @@ subroutine crest_oloop_pr_progress(env,total,current) !* A subroutine to print and track progress of !* concurrent geometry optimizations !********************************************* - use crest_parameters, only:wp,stdout + use crest_parameters,only:wp,stdout use crest_data - use iomod, only: to_str + use iomod,only:to_str implicit none type(systemdata),intent(inout) :: env integer,intent(in) :: total,current @@ -378,24 +440,23 @@ subroutine crest_oloop_pr_progress(env,total,current) call printprogbar(percent) end if increment = 10.0_wp - if(total > 1000) increment = 7.5_wp - if(total > 5000) increment = 5.0_wp - if(total >10000) increment = 2.5_wp - if(total >20000) increment = 1.0_wp - + if (total > 1000) increment = 7.5_wp + if (total > 5000) increment = 5.0_wp + if (total > 10000) increment = 2.5_wp + if (total > 20000) increment = 1.0_wp - else if (current <= total .and. current > 0) then !> the regular printout case + else if (current <= total.and.current > 0) then !> the regular printout case if (env%niceprint) then call printprogbar(percent) else if (.not.env%legacy) then - if(percent >= progressbarrier)then - write(atmp,'(f5.1)') percent - write (stdout,'(1x,a)',advance='no') '|>'//trim(adjustl(atmp))//'%' - progressbarrier = progressbarrier + increment - progressbarrier = min(progressbarrier, 100.0_wp) - flush(stdout) - endif + if (percent >= progressbarrier) then + write (atmp,'(f5.1)') percent + write (stdout,'(1x,a)',advance='no') '|>'//trim(adjustl(atmp))//'%' + progressbarrier = progressbarrier+increment + progressbarrier = min(progressbarrier,100.0_wp) + flush (stdout) + end if else write (stdout,'(1x,i0)',advance='no') current flush (stdout) @@ -429,7 +490,7 @@ subroutine crest_search_multimd(env,mol,mddats,nsim) use dynamics_module use iomod,only:makedir,directory_exist,remove use omp_lib - use crest_restartlog, only: trackrestart,restart_write_dummy + use crest_restartlog,only:trackrestart,restart_write_dummy implicit none type(systemdata),intent(inout) :: env type(mddata) :: mddats(nsim) @@ -450,10 +511,10 @@ subroutine crest_search_multimd(env,mol,mddats,nsim) type(timer) :: profiler !===========================================================! !>--- decide wether to skip this call - if(trackrestart(env))then - call restart_write_dummy('crest_dynamics.trj') - return - endif + if (trackrestart(env)) then + call restart_write_dummy('crest_dynamics.trj') + return + end if !>--- check if we have any MD & calculation settings allocated if (.not.env%mddat%requested) then @@ -465,7 +526,7 @@ subroutine crest_search_multimd(env,mol,mddats,nsim) end if !>--- prepare calculation containers for parallelization (one per thread) - call new_ompautoset(env,'auto_nested',nsim,T,Tn) + call new_ompautoset(env,'auto_nested',nsim,T,Tn) nested = env%omp_allow_nested allocate (calculations(T),source=env%calc) @@ -507,7 +568,7 @@ subroutine crest_search_multimd(env,mol,mddats,nsim) vz = i !>--- OpenMP nested region threads - if(nested) call ompmklset(Tn) + if (nested) call ompmklset(Tn) !!$omp task firstprivate( vz ) private( job,thread_id,io,ex ) call initsignal() @@ -629,11 +690,11 @@ subroutine crest_search_multimd_init(env,mol,mddat,nsim) call move_alloc(calc%calcs(1)%wbo,shk%wbo) end if - if(calc%nfreeze > 0)then + if (calc%nfreeze > 0) then shk%freezeptr => calc%freezelist else - nullify(shk%freezeptr) - endif + nullify (shk%freezeptr) + end if shk%shake_mode = env%shake mddat%shk = shk @@ -708,13 +769,13 @@ subroutine crest_search_multimd_init2(env,mddats,nsim) allocate (mddats(i)%mtd(1),source=mtds(i)) allocate (mddats(i)%cvtype(1),source=cv_rmsd) !> if necessary exclude atoms from RMSD bias - if(sum(env%includeRMSD) /= env%ref%nat)then - if(.not.allocated(mddats(i)%mtd(1)%atinclude)) & - & allocate (mddats(i)%mtd(1)%atinclude( env%ref%nat ), source=.true. ) - do j=1,env%ref%nat - if(env%includeRMSD(j) .ne.1) mddats(i)%mtd(1)%atinclude(j) = .false. - enddo - endif + if (sum(env%includeRMSD) /= env%ref%nat) then + if (.not.allocated(mddats(i)%mtd(1)%atinclude)) & + & allocate (mddats(i)%mtd(1)%atinclude(env%ref%nat),source=.true.) + do j = 1,env%ref%nat + if (env%includeRMSD(j) .ne. 1) mddats(i)%mtd(1)%atinclude(j) = .false. + end do + end if end if end do if (allocated(mtds)) deallocate (mtds) @@ -736,7 +797,7 @@ subroutine crest_search_multimd2(env,mols,mddats,nsim) use shake_module use iomod,only:makedir,directory_exist,remove use omp_lib - use crest_restartlog, only: trackrestart,restart_write_dummy + use crest_restartlog,only:trackrestart,restart_write_dummy implicit none !> INPUT type(systemdata),intent(inout) :: env @@ -756,10 +817,10 @@ subroutine crest_search_multimd2(env,mols,mddats,nsim) type(timer) :: profiler !===========================================================! !>--- decide wether to skip this call - if(trackrestart(env))then - call restart_write_dummy('crest_dynamics.trj') - return - endif + if (trackrestart(env)) then + call restart_write_dummy('crest_dynamics.trj') + return + end if !>--- check if we have any MD & calculation settings allocated if (.not.env%mddat%requested) then @@ -805,7 +866,7 @@ subroutine crest_search_multimd2(env,mols,mddats,nsim) vz = i !>--- OpenMP nested region threads - if(nested) call ompmklset(Tn) + if (nested) call ompmklset(Tn) !$omp task firstprivate( vz ) private( job,thread_id,io,ex ) call initsignal() @@ -914,9 +975,9 @@ subroutine parallel_md_block_printout(MD,vz) write (stdout,'(2x,"| SHAKE algorithm :",a5," (H only) |")') to_str(MD%shake) end if end if - if(allocated(MD%active_potentials))then - write (stdout,'(2x,"| active potentials :",i4," potential |")') size(MD%active_potentials,1) - endif + if (allocated(MD%active_potentials)) then + write (stdout,'(2x,"| active potentials :",i4," potential |")') size(MD%active_potentials,1) + end if if (MD%simtype == type_mtd) then if (MD%cvtype(1) == cv_rmsd) then write (stdout,'(2x,"| dump interval(Vbias) :",f8.2," ps |")') & diff --git a/src/algos/protonate.f90 b/src/algos/protonate.f90 index 26cc0e4b..e7f4e945 100644 --- a/src/algos/protonate.f90 +++ b/src/algos/protonate.f90 @@ -35,6 +35,7 @@ subroutine crest_new_protonate(env,tim) use crest_calculator use strucrd use optimize_module + use parallel_interface implicit none type(systemdata),intent(inout) :: env type(timer),intent(inout) :: tim diff --git a/src/algos/refine.f90 b/src/algos/refine.f90 index c335426d..71a96b28 100644 --- a/src/algos/refine.f90 +++ b/src/algos/refine.f90 @@ -31,6 +31,7 @@ subroutine crest_refine(env,input,output) use crest_calculator use strucrd use cregen_interface + use parallel_interface implicit none type(systemdata),intent(inout) :: env character(len=*),intent(in) :: input diff --git a/src/algos/search_1.f90 b/src/algos/search_1.f90 index c6818082..bd2e8b0e 100644 --- a/src/algos/search_1.f90 +++ b/src/algos/search_1.f90 @@ -27,6 +27,7 @@ subroutine crest_search_1(env,tim) use strucrd use dynamics_module use shake_module + use parallel_interface implicit none type(systemdata),intent(inout) :: env type(timer),intent(inout) :: tim diff --git a/src/algos/search_conformers.f90 b/src/algos/search_conformers.f90 index 9dc00042..dc8d58c1 100644 --- a/src/algos/search_conformers.f90 +++ b/src/algos/search_conformers.f90 @@ -294,6 +294,7 @@ subroutine crest_multilevel_oloop(env,ensnam,multilevel_in) use optimize_module use utilities use crest_restartlog + use parallel_interface implicit none type(systemdata) :: env character(len=*),intent(in) :: ensnam diff --git a/src/algos/search_mecp.f90 b/src/algos/search_mecp.f90 index 5270d5d2..fdf4a91a 100644 --- a/src/algos/search_mecp.f90 +++ b/src/algos/search_mecp.f90 @@ -7,6 +7,7 @@ subroutine crest_search_mecp(env,tim) use dynamics_module use shake_module use cregen_interface + use parallel_interface implicit none type(systemdata),intent(inout) :: env type(timer),intent(inout) :: tim diff --git a/src/qmhelpers/local.f90 b/src/qmhelpers/local.f90 index d4d8cce9..eaca0953 100644 --- a/src/qmhelpers/local.f90 +++ b/src/qmhelpers/local.f90 @@ -99,7 +99,7 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z, & character(len=80) :: atmp character(len=5) :: lmostring(4) data lmostring/'sigma','LP ','pi ','delpi'/ - !data lmostring/'σ ','LP ','π ','δπ '/ + !data lmostring/'σ ','LP ','π ','del.π'/ logical l1,l2,l3,flip integer LWORK,IWORK,LIWORK,INFO integer :: iscreen,icoord,ilmoi,icent ! file handles @@ -276,7 +276,7 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z, & allocate (rklmo(5,2*n)) !if (pr_local) write (*,*) 'lmo centers(Z=2) and atoms on file ' - if (pr_local) write (*,*) 'LMO Fii/eV ncent charge center contributions...' + if (pr_local) write (*,*) ' LMO type Fii/eV ncent charge center contributions...' ! if (pr_local) open (newunit=iscreen,file='xtbscreen.xyz') allocate (tmpq(nat,n)) tmpq(1:nat,1:n) = qmo(1:nat,1:n) @@ -322,7 +322,7 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z, & & imem(1),imem(2),xcen(i),.true.,pithr,jdum) end if if (pr_local) then - write (*,'(i5,1x,a5,1x,2f7.2,3f10.5,12(i5,2x,a2,'':'',f6.2))') & + write (*,'(i5,1x,a,1x,2f7.2,3f10.5,12(i5,2x,a2,'':'',f6.2))') & & i,lmostring(jdum),autoev*f(i),xcen(i),ecent(i,1:3), & & (imem(j),i2e(at(imem(j))),qmo(j,i),j=1,idum) end if @@ -427,7 +427,7 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z, & end do end do - if (pr_local) write (*,*) 'thr ',pithr,'# pi deloc LMO',npi + if (pr_local) write (*,'(1x,a,f6.2,a,i0)') 'ncent threshold ',pithr,' --> # pi deloc LMO ',npi allocate (wbo(nat,nat)) wbo = 0.0d0 @@ -512,12 +512,13 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z, & if(rklmo(5,i) > 1) nprot=nprot+1 enddo if(nprot > 0)then - allocate(protxyz(3,nprot), source=0.0_wp) + allocate(protxyz(4,nprot), source=0.0_wp) m=0 do i=1,k - if(rklmo(5,i) > 1)then + if(nint(rklmo(5,i)) > 1)then m=m+1 protxyz(1:3,m) = rklmo(1:3,i) + protxyz(4,m) = rklmo(5,i) endif enddo endif From c8fdebb0c5eee2fc028e0f7a3e0d37ec7482312b Mon Sep 17 00:00:00 2001 From: Philipp Pracht Date: Fri, 19 Jul 2024 02:01:58 +0200 Subject: [PATCH 3/5] updated gfn0 subproject to handle monoatomic inputs --- subprojects/gfn0 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/subprojects/gfn0 b/subprojects/gfn0 index b0d68ec6..584cec4b 160000 --- a/subprojects/gfn0 +++ b/subprojects/gfn0 @@ -1 +1 @@ -Subproject commit b0d68ec6b44a176db6c3684d7ccc9776e9a50394 +Subproject commit 584cec4b47da23bf3634ef0dd798a1639fcc5e47 From b231aa032acfe668ba55c9eb33b8cb2644160fd1 Mon Sep 17 00:00:00 2001 From: Philipp Pracht Date: Fri, 19 Jul 2024 03:36:59 +0200 Subject: [PATCH 4/5] restore diatomic safety termination for new implementations Signed-off-by: Philipp Pracht --- src/algos/protonate.f90 | 10 +++++ src/algos/search_1.f90 | 6 +++ src/algos/search_conformers.f90 | 6 +++ src/algos/search_entropy.f90 | 7 +++- src/algos/search_newnci.f90 | 6 +++ src/algos/setuptest.f90 | 17 +++++--- src/classes.f90 | 3 +- src/legacy_algos/confscript2_misc.f90 | 4 +- src/legacy_wrappers.f90 | 60 +++++++++++++++++++++------ 9 files changed, 97 insertions(+), 22 deletions(-) diff --git a/src/algos/protonate.f90 b/src/algos/protonate.f90 index e7f4e945..40bb81ca 100644 --- a/src/algos/protonate.f90 +++ b/src/algos/protonate.f90 @@ -36,6 +36,7 @@ subroutine crest_new_protonate(env,tim) use strucrd use optimize_module use parallel_interface + use cregen_interface implicit none type(systemdata),intent(inout) :: env type(timer),intent(inout) :: tim @@ -159,6 +160,15 @@ subroutine crest_new_protonate(env,tim) call rename(ensemblefile,'protonate_1.xyz') + + + +!>--- sorting + write(stdout,'(a)') '> Sorting structures by energy ...' + call newcregen(env,6,'protonate_1.xyz') + write(stdout,'(a)') '> sorted file was renamed to protonated.xyz' + call rename('protonate_1.xyz.sorted','protonated.xyz') + !========================================================================================! return end subroutine crest_new_protonate diff --git a/src/algos/search_1.f90 b/src/algos/search_1.f90 index bd2e8b0e..92899b02 100644 --- a/src/algos/search_1.f90 +++ b/src/algos/search_1.f90 @@ -60,6 +60,12 @@ subroutine crest_search_1(env,tim) call mol%append(stdout) write (stdout,*) +!>--- saftey termination + if(mol%nat .le. 2)then + call catchdiatomic(env) + return + endif + !===========================================================! !>--- Dynamics diff --git a/src/algos/search_conformers.f90 b/src/algos/search_conformers.f90 index dc8d58c1..d4a0b387 100644 --- a/src/algos/search_conformers.f90 +++ b/src/algos/search_conformers.f90 @@ -74,6 +74,12 @@ subroutine crest_search_imtdgc(env,tim) call mol%append(stdout) write (stdout,*) +!>--- saftey termination + if(mol%nat .le. 2)then + call catchdiatomic(env) + return + endif + !>--- sets the MD length according to a flexibility measure call md_length_setup(env) !>--- create the MD calculator saved to env diff --git a/src/algos/search_entropy.f90 b/src/algos/search_entropy.f90 index 028b81bc..1a8e25a5 100644 --- a/src/algos/search_entropy.f90 +++ b/src/algos/search_entropy.f90 @@ -22,7 +22,6 @@ subroutine crest_search_entropy(env,tim) !* This is the re-implementation of CREST's sMTD-iMTD workflow !* from https://doi.org/10.1039/d1sc00621e !* with calculation of conformational entropy -!* This is a TODO !******************************************************************* use crest_parameters,only:wp,stdout use crest_data @@ -81,6 +80,12 @@ subroutine crest_search_entropy(env,tim) call mol%append(stdout) write (stdout,*) +!>--- saftey termination + if(mol%nat .le. 2)then + call catchdiatomic(env) + return + endif + !>--- sets the MD length according to a flexibility measure call md_length_setup(env) !>--- create the MD calculator saved to env diff --git a/src/algos/search_newnci.f90 b/src/algos/search_newnci.f90 index bc6b88c5..f8b1d45e 100644 --- a/src/algos/search_newnci.f90 +++ b/src/algos/search_newnci.f90 @@ -71,6 +71,12 @@ subroutine crest_search_newnci(env,tim) call mol%append(stdout) write (stdout,*) +!>--- saftey termination + if(mol%nat .le. 2)then + call catchdiatomic(env) + return + endif + !>--- sets the MD length according to a flexibility measure call md_length_setup(env) !>--- create the MD calculator saved to env diff --git a/src/algos/setuptest.f90 b/src/algos/setuptest.f90 index 78132d9e..baf35ecb 100644 --- a/src/algos/setuptest.f90 +++ b/src/algos/setuptest.f90 @@ -273,6 +273,7 @@ subroutine trialOPT_calculator(env) use crest_calculator use optimize_module use strucrd + use iomod implicit none !> INPUT type(systemdata),intent(inout) :: env @@ -301,11 +302,19 @@ subroutine trialOPT_calculator(env) !>--- perform geometry optimization pr = .false. !> stdout printout wr = .true. !> write crestopt.log + if(wr)then + call remove('crestopt.log') + endif call optimize_geometry(mol,molopt,tmpcalc,energy,grd,pr,wr,io) !>--- check success success = (io == 0) call trialOPT_warning(env,molopt,success) +!>--- if the checks were successfull, env%ref is overwritten + env%ref%nat = molopt%nat + env%ref%at = molopt%at + env%ref%xyz = molopt%xyz + env%ref%etot = energy deallocate(grd) end subroutine trialOPT_calculator @@ -336,7 +345,7 @@ subroutine trialOPT_warning(env,mol,success) if (.not.success) then write (stdout,*) write (stdout,*) ' Initial geometry optimization failed!' - write (stdout,*) ' Please check your input.' + write (stdout,*) ' Please check your input and, if present, crestopt.log.' error stop end if write (stdout,*) 'Geometry successfully optimized.' @@ -395,11 +404,7 @@ subroutine trialOPT_warning(env,mol,success) end if end if -!>--- If all checks succeded, update reference with optimized geometry - env%ref%nat = mol%nat - env%ref%at = mol%at - env%ref%xyz = mol%xyz - + return end subroutine trialOPT_warning !========================================================================================! diff --git a/src/classes.f90 b/src/classes.f90 index 666d96dc..e2dd4dde 100644 --- a/src/classes.f90 +++ b/src/classes.f90 @@ -272,7 +272,8 @@ module crest_data real(wp),allocatable :: xyz(:,:) integer :: ichrg = 0 integer :: uhf = 0 - integer :: ntopo + integer :: ntopo = 0 + real(wp) :: etot = 0.0_wp integer,allocatable :: topo(:) real(wp),allocatable :: charges(:) real(wp),allocatable :: wbo(:,:) diff --git a/src/legacy_algos/confscript2_misc.f90 b/src/legacy_algos/confscript2_misc.f90 index 3a8bffc0..b92302d5 100644 --- a/src/legacy_algos/confscript2_misc.f90 +++ b/src/legacy_algos/confscript2_misc.f90 @@ -1147,7 +1147,7 @@ end subroutine elowcheck !----------------------------------------------------------------------- ! diatomic "error" catcher ---> we will have no conformers !----------------------------------------------------------------------- -subroutine catchdiatomic(env) +subroutine catchdiatomic_legacy(env) use crest_parameters use crest_data use strucrd @@ -1174,7 +1174,7 @@ subroutine catchdiatomic(env) call mol%deallocate return -end subroutine catchdiatomic +end subroutine catchdiatomic_legacy !-------------------------------------------------------------------------------------! ! entropy ensemble file copy routine diff --git a/src/legacy_wrappers.f90 b/src/legacy_wrappers.f90 index 83d88899..c5825d4c 100644 --- a/src/legacy_wrappers.f90 +++ b/src/legacy_wrappers.f90 @@ -52,11 +52,11 @@ subroutine env2calc(env,calc,molin) cal%rdwbo = .false. cal%rddip = .true. !> except for SP runtype (from command line!) - if( env%crestver == crest_sp )then + if (env%crestver == crest_sp) then cal%rdwbo = .true. cal%rddip = .true. cal%rdqat = .true. - endif + end if !> implicit solvation if (env%gbsa) then @@ -89,11 +89,11 @@ subroutine env2calc(env,calc,molin) end if call cal2%autocomplete(2) - + cal2%refine_lvl = refine%singlepoint call calc%add(cal2) - if(allocated(env%refine_queue)) deallocate(env%refine_queue) - call env%addrefine( refine%singlepoint ) + if (allocated(env%refine_queue)) deallocate (env%refine_queue) + call env%addrefine(refine%singlepoint) end if return @@ -127,11 +127,11 @@ end subroutine env2calc end interface - call env%ref%to(mol) + call env%ref%to(mol) call env2calc(env,env%calc,mol) - ! env%calc = calc + ! env%calc = calc end subroutine env2calc_setup !================================================================================! @@ -330,7 +330,7 @@ subroutine trialOPT(env) !* saved to env%ref and checks for changes in the topology !********************************************************** use crest_data - use crest_parameters, only: stdout + use crest_parameters,only:stdout implicit none !> INPUT type(systemdata) :: env @@ -341,15 +341,21 @@ subroutine trialOPT(env) call trialOPT_calculator(env) end if - if(env%crestver == crest_trialopt)then + if (env%crestver == crest_trialopt) then !>-- if we reach this point in the standalone trialopt the geometry is ok! - write(stdout,*) - stop 'Geometry ok!' - endif + write (stdout,*) + stop 'Geometry ok!' + end if end subroutine trialOPT !================================================================================! + subroutine protonate(env,tim) +!***************************************************** +!* subroutine protonate +!* driver for the automated protonation site search +!* originally published in JCC +!***************************************************** use iso_fortran_env,only:wp => real64 use crest_data implicit none @@ -362,4 +368,34 @@ subroutine protonate(env,tim) end if end subroutine protonate +!========================================================================================! +subroutine catchdiatomic(env) +!**************************************** +!* subroutine catchdiatomic +!* if we only have one or two atoms just +!* write the "optimized" structure +!**************************************** + use crest_data + use crest_parameters + use strucrd + use iomod,only:copy + use cregen_interface + implicit none + type(systemdata) :: env + integer :: ich + type(coord) :: mol + if (env%legacy) then + call catchdiatomic_legacy(env) + else + call env%ref%to(mol) + open (file=conformerfile,newunit=ich) + mol%xyz = mol%xyz*bohr !to ang + call wrxyz(ich,mol%nat,mol%at,mol%xyz,env%ref%etot) + close (ich) + call copy('xtbopt.xyz',conformerfile) + call copy(conformerfile,'crest_rotamers.xyz') + call copy(conformerfile,'crest_best.xyz') + end if + call newcregen(env,6,'crest_rotamers.xyz') +end subroutine catchdiatomic From 31fede016b43c55de6196bdabec554a30e8af691 Mon Sep 17 00:00:00 2001 From: Philipp Pracht Date: Fri, 19 Jul 2024 03:58:58 +0200 Subject: [PATCH 5/5] default back to old protonation search for now --- src/legacy_wrappers.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/legacy_wrappers.f90 b/src/legacy_wrappers.f90 index c5825d4c..b0ddd427 100644 --- a/src/legacy_wrappers.f90 +++ b/src/legacy_wrappers.f90 @@ -361,11 +361,11 @@ subroutine protonate(env,tim) implicit none type(systemdata) :: env type(timer) :: tim - if (env%legacy) then + !if (env%legacy) then call protonate_legacy(env,tim) - else - call crest_new_protonate(env,tim) - end if + !else + ! call crest_new_protonate(env,tim) + !end if end subroutine protonate !========================================================================================!