diff --git a/src/algos/parallel.f90 b/src/algos/parallel.f90 index 13fd6f10..71afe979 100644 --- a/src/algos/parallel.f90 +++ b/src/algos/parallel.f90 @@ -225,16 +225,18 @@ end subroutine crest_sploop !========================================================================================! !========================================================================================! 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 +!* WARNING: the ensemble file will NOT be in the same order +!* as the input xyz array. However, the overwritten xyz will be! !* customcalc - customized (optional) calculation level data !* !* IMPORTANT: xyz should be in Bohr(!) for this routine -!*************************************************************** +!****************************************************************************** use crest_parameters,only:wp,stdout,sep use crest_calculator use omp_lib @@ -369,6 +371,11 @@ subroutine crest_oloop(env,nat,nall,at,xyz,eread,dump,customcalc) end if eread(zcopy) = energy xyz(:,:,zcopy) = molsnew(job)%xyz(:,:) + else if(io==calculations(job)%maxcycle .and. calculations(job)%anopt) then + !>--- allow partial optimization? + c = c+1 + eread(zcopy) = energy + xyz(:,:,zcopy) = molsnew(job)%xyz(:,:) else eread(zcopy) = 1.0_wp end if diff --git a/src/algos/protonate.f90 b/src/algos/protonate.f90 index 40bb81ca..62c564c2 100644 --- a/src/algos/protonate.f90 +++ b/src/algos/protonate.f90 @@ -37,6 +37,7 @@ subroutine crest_new_protonate(env,tim) use optimize_module use parallel_interface use cregen_interface + use iomod implicit none type(systemdata),intent(inout) :: env type(timer),intent(inout) :: tim @@ -51,12 +52,14 @@ subroutine crest_new_protonate(env,tim) real(wp) :: energy,dip real(wp),allocatable :: grad(:,:) type(calcdata),allocatable :: tmpcalc + type(calcdata),allocatable :: tmpcalc_ff + type(calculation_settings) :: tmpset real(wp),allocatable :: protxyz(:,:) - integer :: natp + integer :: natp,pstep,npnew integer,allocatable :: atp(:) - real(wp),allocatable :: xyzp(:,:,:) + real(wp),allocatable :: xyzp(:,:,:) real(wp),allocatable :: ep(:) - character(len=*),parameter :: partial = '∂E/∂' + logical,allocatable :: atlist(:) !========================================================================================! write (stdout,*) !call system('figlet singlepoint') @@ -68,10 +71,11 @@ subroutine crest_new_protonate(env,tim) write (stdout,*) "|_| " write (stdout,*) "-----------------------------------------------" write (stdout,*) " automated protonation site screening script " + write (stdout,*) " revised version (c) P.Pracht 2024 " write (stdout,*) 'Cite as:' - write (stdout,*) 'P.Pracht, C.A.Bauer, S.Grimme' - write (stdout,*) 'JCC, 2017, 38, 2618–2631.' - write (stdout,*) + 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) @@ -84,6 +88,8 @@ subroutine crest_new_protonate(env,tim) write (stdout,*) !========================================================================================! !>--- The first step is to perfom a LMO calculation to identify suitable protonation sites +!>--- in this version of the program this step is always done with GFN0-xTB. + call tim%start(14,'LMO center calculation') write (stdout,'(a)') repeat('-',80) write (stdout,'(a)') @@ -91,11 +97,13 @@ subroutine crest_new_protonate(env,tim) flush (stdout) allocate (grad(3,mol%nat),source=0.0_wp) +!>--- GFN0 job adapted from global settings allocate (tmpcalc) call env2calc(env,tmpcalc,mol) tmpcalc%calcs(1)%id = jobtype%gfn0 tmpcalc%calcs(1)%rdwbo = .true. tmpcalc%calcs(1)%getlmocent = .true. + tmpcalc%calcs(1)%chrg = env%chrg tmpcalc%ncalculations = 1 write (stdout,'(a)') 'done.' !>--- and then start it @@ -119,15 +127,25 @@ subroutine crest_new_protonate(env,tim) !>--- 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' + 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 + write (stdout,'(1x,a5,1x,a,5x,a)') 'LMO','type','center(xyz/Ang)' + do i = 1,np + select case (nint(protxyz(4,i))) + case (1) + write (stdout,'(1x,i5,1x,a,3F12.5)') i,'LP ',protxyz(1:3,i)*autoaa + case (2) + write (stdout,'(1x,i5,1x,a,3F12.5)') i,'π ',protxyz(1:3,i)*autoaa + case (3) + write (stdout,'(1x,i5,1x,a,3F12.5)') i,'del.π',protxyz(1:3,i)*autoaa + case default + write (stdout,'(1x,i5,1x,a,3F12.5)') i,'??? ',protxyz(1:3,i)*autoaa + end select + end do 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,*) ' Confirm whether you expect π- or LP-centers for your molecule!' write (stdout,*) return end if @@ -136,39 +154,156 @@ subroutine crest_new_protonate(env,tim) !========================================================================================! !>--- If we reached this point, we have candidate positions for our protonation! - write(stdout,'(a)',advance='yes') '> Generating candidate structures ... ' - flush(stdout) + pstep = 0 + 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,*) + allocate (atp(natp),source=0) + allocate (xyzp(3,natp,np),ep(np),source=0.0_wp) -!>--- Enforce further constraints, conditions, etc. -! (TODO) + call protonation_candidates(env,mol,natp,np,protxyz,atp,xyzp,npnew) +!>--- NOTE: after this the global charge (env%chrg) and all charges saved in the calc levels have been increased + + if (npnew < 1) then + write (stdout,*) + write (stdout,'(a)') '> WARNING: No remaining protonation sites after applying user defined conditions!' + write (stdout,'(a)') '> Modify the search criteria and check your input structure for sanity.' + return + end if + + write (atmp,'(a,i0,a)') 'protonate_',pstep,'.xyz' + write (stdout,'(a,a,a,i0,a)') '> Write ',trim(atmp),' with ',npnew,' candidates ... ' -!>--- 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') + call wrensemble('protonate_0.xyz',natp,npnew,atp,xyzp(:,:,1:npnew)*autoaa,ep(1:npnew)) + write (stdout,'(a)') '> done.' + write (stdout,*) +!========================================================================================! +!>--- Optimize candidates, optional FF pre-step + if (env%protb%ffopt) then + call smallhead('Protomer Ensemble FF Pre-Optimization') + !>--- set up a temporary calculation object + write (stdout,'(a)') '> LMO centers can be very close to atoms, leading to initial' + write (stdout,'(a)') ' extremely high-energy candidates which may undergo unwanted' + write (stdout,'(a)') ' chemical changes in optimizations. Classical force-fields' + write (stdout,'(a)') ' with defined bond-topology can help circumvent this issue.' + write (stdout,'(a)') '> Setting up force-field structure pre-optimization ...' + allocate (tmpcalc_ff) + tmpcalc_ff%optnewinit = .true. + env%calc%optnewinit = .true. + tmpcalc_ff%optnewinit=.true. + call tmpset%create('gfnff') + tmpset%chrg = env%chrg + call tmpcalc_ff%add(tmpset) + tmpcalc_ff%maxcycle = 10000 + call tmpcalc_ff%info(stdout) + !>--- Run optimizations + call tim%start(15,'Ensemble optimization (FF)') + call print_opt_data(tmpcalc_ff,stdout) + write (stdout,'(a,i0,a)') '> ',npnew,' structures to optimize ...' + call crest_oloop(env,natp,npnew,atp,xyzp(:,:,1:npnew),ep(1:npnew),.false.,tmpcalc_ff) + call tim%stop(15) + deallocate (tmpcalc_ff) -!>--- 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') + pstep = pstep+1 + write (atmp,'(a,i0,a)') 'protonate_',pstep,'.xyz' + write (stdout,'(a,a,a)') '> Write ',trim(atmp),' with optimized structures ... ' + call wrensemble(trim(atmp),natp,npnew,atp,xyzp(:,:,1:npnew)*autoaa,ep(1:npnew)) + + !>--- sorting + write (stdout,'(a)') '> Sorting structures by energy to remove failed opts. ...' + env%ewin = 5000.0_wp + call newcregen(env,7,trim(atmp)) + call rename(trim(atmp)//'.sorted',trim(atmp)) + write (stdout,'(a)') '> WARNING: These are force-field energies which are ' + write (stdout,'(a)') ' NOT(!) accurate for bond formation and breaking!' + write (stdout,*) + + !>--- re-read sorted ensemble + deallocate (xyzp,atp) !> clear this space to re-use it + call rdensemble(trim(atmp),natp,npnew,atp,xyzp) + xyzp = xyzp*aatoau !> don't forget to restore BOHR + end if +!========================================================================================! +!>--- H-position-only optimization (only makes sense after FF preoptimization) + if (env%protb%hnewopt.and.env%protb%ffopt) then + call smallhead('Protomer Ensemble Frozen-Atom Optimization') + !>--- create temporary calculation from our intended level of theory + write (stdout,'(a)') '> Setting up frozen structure optimization ...' + allocate (tmpcalc) + call env2calc(env,tmpcalc,mol) + !>--- freeze all atoms, EXCEPT the new one + allocate (atlist(natp),source=.true.) + atlist(natp) = .false. !> the new one is always last + do i = 1,natp + if (atp(i) == 1) then + atlist(i) = .false. !> additionally un-freeze all H's (this seems to be beneficial) + end if + end do + tmpcalc%nfreeze = count(atlist) + write (stdout,'(a,i0,a)') '> ',tmpcalc%nfreeze,' frozen atoms. All H non-frozen.' + call move_alloc(atlist,tmpcalc%freezelist) + call tmpcalc%info(stdout) + + !>--- run opt + call tim%start(16,'Ensemble optimization (frozen)') + write (stdout,'(a,i0,a)') '> ',npnew,' structures to optimize ...' + call crest_oloop(env,natp,npnew,atp,xyzp(:,:,1:npnew),ep(1:npnew),.false.,tmpcalc) + call tim%stop(16) + + pstep = pstep+1 + write (atmp,'(a,i0,a)') 'protonate_',pstep,'.xyz' + write (stdout,'(a,a,a)') '> Write ',trim(atmp),' with optimized structures ... ' + call wrensemble(trim(atmp),natp,npnew,atp,xyzp(:,:,1:npnew)*autoaa,ep(1:npnew)) + + !>--- sorting + write (stdout,'(a)') '> Sorting structures by energy to remove failed opts. ...' + env%ewin = env%protb%ewin*10.0_wp !> large energy threshold + call newcregen(env,7,trim(atmp)) + call rename(trim(atmp)//'.sorted',trim(atmp)) + write (stdout,*) + + !>--- re-read sorted ensemble + deallocate (xyzp,atp) !> clear this space to re-use it + call rdensemble(trim(atmp),natp,npnew,atp,xyzp) + xyzp = xyzp*aatoau !> don't forget to restore BOHR + + deallocate (tmpcalc) + end if + +!========================================================================================! +!>--- Optimize with global settings + if (env%protb%finalopt.and.env%protb%ffopt) then + call smallhead('Final Protomer Ensemble Optimization') + allocate(tmpcalc) + call env2calc(env,tmpcalc,mol) + call tmpcalc%info(stdout) + !tmpcalc%maxcycle=5 + !tmpcalc%anopt=.true. + call tim%start(20,'Ensemble optimization') + call print_opt_data(env%calc,stdout) + write (stdout,'(a,i0,a)') '> ',npnew,' structures to optimize ...' + call crest_oloop(env,natp,npnew,atp,xyzp(:,:,1:npnew),ep,.false.,tmpcalc) + call tim%stop(20) + + pstep = pstep+1 + write (atmp,'(a,i0,a)') 'protonate_',pstep,'.xyz' + write (stdout,'(a,a,a)') '> Write ',trim(atmp),' with optimized structures ... ' + call wrensemble(trim(atmp),natp,npnew,atp,xyzp(:,:,1:npnew)*autoaa,ep(1:npnew)) + +!>--- sorting + write (stdout,'(a)') '> Sorting structures by energy ...' + env%ewin = env%protb%ewin + call newcregen(env,7,trim(atmp)) + call rename(trim(atmp)//'.sorted','protonated.xyz') + else + call rename(trim(atmp),'protonated.xyz') + end if + write (stdout,'(a)') '> sorted file was renamed to protonated.xyz' !========================================================================================! return end subroutine crest_new_protonate @@ -176,7 +311,7 @@ end subroutine crest_new_protonate !========================================================================================! !========================================================================================! -subroutine protonation_candidates(env,mol,natp,np,protxyz,at,xyz) +subroutine protonation_candidates(env,mol,natp,np,protxyz,at,xyz,npnew) !******************************************************** !* generate protonation/ionization candidate structures !* The outputs are at and xyz, the latter being in Bohr @@ -190,30 +325,60 @@ subroutine protonation_candidates(env,mol,natp,np,protxyz,at,xyz) type(coord),intent(in) :: mol integer,intent(in) :: natp integer,intent(in) :: np - real(wp),intent(in) :: protxyz(3,np) + real(wp),intent(in) :: protxyz(4,np) !> OUTPUT integer,intent(out) :: at(natp) real(wp),intent(out) :: xyz(3,natp,np) + integer,intent(out) :: npnew !> LOCAL - integer :: i,j,k,l + integer :: i,j,k,l,ii + integer :: ati,ichrg,ctype 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() + if (env%protb%swelem) then +!>--- User-defined monoatomic ion + ichrg = env%protb%swchrg + ati = env%protb%swat + else +!>--- DEFAULT: H⁺ + ichrg = 1 + ati = 1 + endif + write (stdout,'(a,i0)') '> Increasing the molecular charge by ',ichrg + call env%calc%increase_charge(ichrg) + env%chrg = env%chrg + ichrg + +!>--- Check if we have some other conditions + if(any(.not.env%protb%active_lmo(:)))then + write(stdout,'(a)',advance='no') '> User-defined: IGNORING ' + if(.not.env%protb%active_lmo(1)) write(stdout,'(a)',advance='no') 'π ' + if(.not.env%protb%active_lmo(2)) write(stdout,'(a)',advance='no') 'LP ' + if(.not.env%protb%active_lmo(3)) write(stdout,'(a)',advance='no') 'deloc.π ' + write(stdout,'(a)') 'LMOs ...' + endif !>--- Populate + npnew = 0 + ii = 0 do i = 1,np +!>--- Enforce further constraints, conditions, etc. + ctype = nint(protxyz(4,i)) + if(.not.env%protb%active_lmo(ctype)) cycle !> skip unselected LMO types + + ii= ii + 1 !> counter of actually created structures do j = 1,mol%nat - xyz(1:3,j,i) = mol%xyz(1:3,j) + xyz(1:3,j,ii) = mol%xyz(1:3,j) at(j) = mol%at(j) end do - xyz(1:3,natp,i) = protxyz(1:3,i) - at(natp) = 1 + xyz(1:3,natp,ii) = protxyz(1:3,i) + at(natp) = ati end do + npnew = ii + end subroutine protonation_candidates !========================================================================================! diff --git a/src/calculator/calc_type.f90 b/src/calculator/calc_type.f90 index 8e48dff4..d87074e7 100644 --- a/src/calculator/calc_type.f90 +++ b/src/calculator/calc_type.f90 @@ -68,10 +68,12 @@ module calc_type !&> !=========================================================================================! -!>--- data object that contains the data for a *SINGLE* calculation + public :: calculation_settings type :: calculation_settings - +!********************************************************************** +!* data object that contains the data for a *SINGLE* calculation level +!********************************************************************** integer :: id = 0 !> calculation type (see "jobtype" parameter above) integer :: prch = stdout !> printout channel logical :: pr = .false. !> allow the calculation to produce printout? Results in a lot I/O @@ -148,7 +150,7 @@ module calc_type character(len=:),allocatable :: tbliteparam !>--- GFN0-xTB data - type(gfn0_data),allocatable :: g0calc + type(gfn0_data),allocatable :: g0calc integer :: nconfig = 0 integer,allocatable :: config(:) real(wp),allocatable :: occ(:) @@ -181,13 +183,16 @@ module calc_type procedure :: norestarts => calculation_settings_norestarts procedure :: dumpdipgrad => calculation_dump_dipgrad end type calculation_settings -!=========================================================================================! !=========================================================================================! -!> data object that collects settings for *ALL* calculations and constraints. +!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< this parameter will decide how to return or add up energies and gradients integer :: refine_stage = 0 !> to allow iterating different refinement stages @@ -223,6 +228,8 @@ module calc_type real(wp) :: etot !>--- optimization settings + logical :: optnewinit = .false. !> ensure fresh calc/param setup at beginning of opt + logical :: anopt = .false. !> allow collecting structures that ran into maxcycle? integer :: optlev = 0 integer :: micro_opt = 20 integer :: maxcycle = 0 @@ -272,6 +279,7 @@ module calc_type procedure :: freezegrad => calculation_freezegrad procedure :: increase_charge => calculation_increase_charge procedure :: decrease_charge => calculation_decrease_charge + procedure :: dealloc_params => calculation_deallocate_params end type calcdata public :: get_dipoles @@ -328,59 +336,25 @@ subroutine calculation_reset(self) return end subroutine calculation_reset - subroutine calculation_settings_deallocate(self) - implicit none - class(calculation_settings) :: self - - if (allocated(self%calcspace)) deallocate (self%calcspace) - if (allocated(self%calcfile)) deallocate (self%calcfile) - if (allocated(self%gradfile)) deallocate (self%gradfile) - if (allocated(self%path)) deallocate (self%path) - if (allocated(self%other)) deallocate (self%other) - if (allocated(self%binary)) deallocate (self%binary) - if (allocated(self%systemcall)) deallocate (self%systemcall) - if (allocated(self%description)) deallocate (self%description) - if (allocated(self%wbo)) deallocate (self%wbo) - if (allocated(self%dipgrad)) deallocate (self%dipgrad) - if (allocated(self%gradkey)) deallocate (self%gradkey) - if (allocated(self%efile)) deallocate (self%efile) - if (allocated(self%solvmodel)) deallocate (self%solvmodel) - if (allocated(self%solvent)) deallocate (self%solvent) - if (allocated(self%tblite)) deallocate (self%tblite) - if (allocated(self%g0calc)) deallocate (self%g0calc) - if (allocated(self%ff_dat)) deallocate (self%ff_dat) - if (allocated(self%xhcff)) deallocate (self%xhcff) - - self%id = 0 - self%prch = stdout - self%chrg = 0 - self%uhf = 0 - self%epot = 0.0_wp - self%efix = 0.0_wp - self%etot = 0.0_wp - - self%rdwbo = .false. - self%rddip = .false. - self%dipole = 0.0_wp - self%rddipgrad = .false. - self%gradtype = 0 - self%gradfmt = 0 - - self%tblitelvl = 2 - self%etemp = 300.0_wp - self%accuracy = 1.0_wp - self%apiclean = .false. - self%maxscc = 500 - self%saveint = .false. - - self%ngrid = 230 - self%extpressure = 0.0_wp - self%proberad = 1.5_wp +!=========================================================================================! - self%ONIOM_highlowroot = 0 - self%ONIOM_id = 0 - return - end subroutine calculation_settings_deallocate + subroutine calculation_deallocate_params(self) +!********************************************** +!* Deallocate parametrization/model setup +!* This only applies (or is useful) for a few +!* selected model like GFN-FF, gfn0 and tblite +!********************************************** + class(calcdata) :: self + integer :: i,j,k + if (self%ncalculations > 0) then + do i = 1,self%ncalculations + if(allocated(self%calcs(i)%tblite)) deallocate(self%calcs(i)%tblite) + if(allocated(self%calcs(i)%g0calc)) deallocate(self%calcs(i)%g0calc) + if(allocated(self%calcs(i)%ff_dat)) deallocate(self%calcs(i)%ff_dat) + if(allocated(self%calcs(i)%xhcff)) deallocate(self%calcs(i)%xhcff) + end do + end if + end subroutine calculation_deallocate_params !=========================================================================================! @@ -660,137 +634,63 @@ end subroutine calculation_decrease_charge !=========================================================================================! - subroutine calculation_settings_addconfig(self,config) - implicit none - class(calculation_settings) :: self - integer,intent(in) :: config(:) - integer :: i,j - integer :: l,lold,lnew,n,nnew - integer,allocatable :: configtmp(:,:) - - l = size(config,1) - if (allocated(self%config)) deallocate (self%config) - allocate (self%config(l)) - self%config = config - end subroutine calculation_settings_addconfig - -!=========================================================================================! - - subroutine calculation_settings_norestarts(self) -!************************************************* -!* remove restart options from this calculation -!************************************************* - implicit none - class(calculation_settings) :: self - self%restart = .false. - if (allocated(self%refgeo)) deallocate (self%refgeo) - if (allocated(self%restartfile)) deallocate (self%restartfile) - end subroutine calculation_settings_norestarts - -!=========================================================================================! - -!>-- check for missing settings in a calculation_settings object - subroutine calculation_settings_autocomplete(self,id) + subroutine calc_set_active(self,ids) implicit none - class(calculation_settings) :: self - integer,intent(in) :: id - integer :: i,j - character(len=50) :: nmbr - - !> add a short description - self%description = trim(jobdescription(self%id+1)) - call calculation_settings_shortflag(self) - - if (.not.allocated(self%calcspace)) then - !> I've decided to perform all calculations in a separate directory to - !> avoid accumulation of files in the main workspace - write (nmbr,'(i0)') id - self%calcspace = 'calculation.level.'//trim(nmbr) - end if - - if (self%pr) then - self%prch = self%prch+id - end if - end subroutine calculation_settings_autocomplete + class(calcdata) :: self + integer,intent(in) :: ids(:) + integer :: i,j,k,l + if (allocated(self%weightbackup)) deallocate (self%weightbackup) +!>--- on-the-fly multiscale definition + allocate (self%weightbackup(self%ncalculations),source=1.0_wp) + do i = 1,self%ncalculations +!>--- save backup weights + self%weightbackup(i) = self%calcs(i)%weight +!>--- set the weight of all unwanted calculations to 0 + if (.not.any(ids(:) .eq. i)) then + self%calcs(i)%weight = 0.0_wp + self%calcs(i)%active = .false. + else +!>--- and all other to 1 + self%calcs(i)%weight = 1.0_wp + self%calcs(i)%active = .true. + end if + end do + end subroutine calc_set_active -!>--- create a short calculation info flag - subroutine calculation_settings_shortflag(self) + subroutine calc_set_active_restore(self) implicit none - class(calculation_settings) :: self - integer :: i,j - - select case( self%id ) - case( jobtype%xtbsys ) - self%shortflag = 'xtb subprocess' - case( jobtype%generic ) - self%shortflag = 'generic subprocess' - case( jobtype%turbomole ) - self%shortflag = 'TURBOMOLE subprocess' - case( jobtype%orca ) - self%shortflag = 'ORCA subprocess' - case( jobtype%terachem ) - self%shortflag = 'TeraChem subprocess' - case( jobtype%tblite ) - select case (self%tblitelvl) - case (xtblvl%gfn2) - self%shortflag = 'GFN2-xTB' - case (xtblvl%gfn1) - self%shortflag = 'GFN1-xTB' - case (xtblvl%ipea1) - self%shortflag = 'IPEA1-xTB' - case (xtblvl%ceh) - self%shortflag = 'CEH' - case (xtblvl%eeq) - self%shortflag = 'EEQ(D4)' - case (xtblvl%param) - self%shortflag = 'parameter file: '//trim(self%tbliteparam) - end select - case( jobtype%gfn0 ) - self%shortflag = 'GFN0-xTB' - case( jobtype%gfn0occ ) - self%shortflag = 'GFN0-xTB*' - case( jobtype%gfnff ) - self%shortflag = 'GFN-FF' - case( jobtype%xhcff ) - self%shortflag = 'XHCFF' - case( jobtype%lj ) - self%shortflag = 'LJ' - case default - self%shortflag = 'undefined' - end select - if(allocated(self%solvmodel).and.allocated(self%solvent))then - self%shortflag = self%shortflag//'/'//trim(self%solvmodel) - self%shortflag = self%shortflag//'('//trim(self%solvent)//')' - endif - end subroutine calculation_settings_shortflag - - + class(calcdata) :: self + integer :: i,j,k,l + if (.not.allocated(self%weightbackup)) return + do i = 1,self%ncalculations +!>--- set all to active and restore saved weights + self%calcs(i)%weight = self%weightbackup(i) + self%calcs(i)%active = .true. + self%eweight(i) = self%weightbackup(i) + end do + deallocate (self%weightbackup) + end subroutine calc_set_active_restore +!==========================================================================================! -!>-- generate a unique print id for the calculation - subroutine calculation_settings_printid(self,thread,id) + subroutine get_dipoles(calc,diptmp) +!********************************************* +!* Collect the x y and z dipole moments for +!* each calculation level in one array diptmp +!********************************************* implicit none - class(calculation_settings) :: self - integer,intent(in) :: thread,id - integer :: i,j,dum - character(len=50) :: nmbr - dum = 100*(thread+1) - dum = dum+id - self%prch = dum - end subroutine calculation_settings_printid + type(calcdata),intent(inout) :: calc + real(wp),allocatable,intent(out) :: diptmp(:,:) + integer :: i,j,k,l,m - subroutine calculation_dump_dipgrad(self,filename) - implicit none - class(calculation_settings) :: self - character(len=*),intent(in) :: filename - integer :: i,j,ich - if (.not.allocated(self%dipgrad)) return - open (newunit=ich,file=filename) - do i = 1,size(self%dipgrad,2) - write (ich,'(3f20.10)') self%dipgrad(1:3,i) + m = calc%ncalculations + allocate (diptmp(3,m),source=0.0_wp) + do i = 1,m + if (calc%calcs(i)%rddip) then + diptmp(1:3,i) = calc%calcs(i)%dipole(1:3) + end if end do - close (ich) - end subroutine calculation_dump_dipgrad + end subroutine get_dipoles !=========================================================================================! @@ -893,106 +793,17 @@ subroutine calculation_ONIOMexpand(self) write (stdout,*) 'done.' end subroutine calculation_ONIOMexpand -!=========================================================================================! - subroutine calculation_settings_info(self,iunit) +!========================================================================================! + + subroutine calculation_info(self,iunit) implicit none - class(calculation_settings) :: self + class(calcdata) :: self integer,intent(in) :: iunit integer :: i,j - character(len=*),parameter :: fmt1 = '(" :",2x,a20," : ",i0)' - character(len=*),parameter :: fmt2 = '(" :",2x,a20," : ",f0.5)' - character(len=*),parameter :: fmt3 = '(" :",2x,a20," : ",a)' - character(len=*),parameter :: fmt4 = '(" :",1x,a)' + character(len=*),parameter :: fmt1 = '(1x,a20," : ",i5)' + character(len=*),parameter :: fmt2 = '(1x,a20," : ",f12.5)' character(len=20) :: atmp - - if (allocated(self%description)) then - write (iunit,'(" :",1x,a)') trim(self%description) - else - write (atmp,*) 'Job type' - write (iunit,fmt1) atmp,self%id - end if - !> more info - if (self%id == jobtype%tblite) then - select case (self%tblitelvl) - case (xtblvl%gfn2) - write (iunit,fmt4) 'GFN2-xTB level' - case (xtblvl%gfn1) - write (iunit,fmt4) 'GFN1-xTB level' - case (xtblvl%ceh) - write (iunit,fmt4) 'Charge Extended Hückel (CEH) model' - end select - end if - if (any((/jobtype%orca,jobtype%xtbsys,jobtype%turbomole, & - & jobtype%generic,jobtype%terachem/) == self%id)) then - write (iunit,'(" :",3x,a,a)') 'selected binary : ',trim(self%binary) - end if - if (self%refine_lvl > 0) then - write (atmp,*) 'refinement stage' - write (iunit,fmt1) atmp,self%refine_lvl - end if - - !> system data - write (atmp,*) 'Molecular charge' - write (iunit,fmt1) atmp,self%chrg - if (self%uhf /= 0) then - write (atmp,*) 'UHF parameter' - write (iunit,fmt1) atmp,self%uhf - end if - - if (allocated(self%solvmodel)) then - write (atmp,*) 'Solvation model' - write (iunit,fmt3) atmp,trim(self%solvmodel) - end if - if (allocated(self%solvent)) then - write (atmp,*) 'Solvent' - write (iunit,fmt3) atmp,trim(self%solvent) - end if - - !> xTB specific parameters - if (any((/jobtype%tblite,jobtype%xtbsys,jobtype%gfn0,jobtype%gfn0occ/) == self%id)) then - write (atmp,*) 'Fermi temperature' - write (iunit,fmt2) atmp,self%etemp - write (atmp,*) 'Accuracy' - write (iunit,fmt2) atmp,self%accuracy - write (atmp,*) 'max SCC cycles' - write (iunit,fmt1) atmp,self%maxscc - end if - - write (atmp,*) 'Reset data?' - if (self%apiclean) write (iunit,fmt3) atmp,'yes' - write (atmp,*) 'Read WBOs?' - if (self%rdwbo) write (iunit,fmt3) atmp,'yes' - write (atmp,*) 'Read dipoles?' - if (self%rddip) write (iunit,fmt3) atmp,'yes' - - if (self%ONIOM_highlowroot /= 0) then - select case (self%ONIOM_highlowroot) - case (1) - write (atmp,*) 'ONIOM frag ("high")' - case (2) - write (atmp,*) 'ONIOM frag ("low")' - case (3) - write (atmp,*) 'ONIOM frag ("root")' - end select - write (iunit,fmt1) trim(atmp),self%ONIOM_id - else - write (atmp,*) 'Weight' - write (iunit,fmt2) atmp,self%weight - end if - - end subroutine calculation_settings_info - -!========================================================================================! - - subroutine calculation_info(self,iunit) - implicit none - class(calcdata) :: self - integer,intent(in) :: iunit - integer :: i,j - character(len=*),parameter :: fmt1 = '(1x,a20," : ",i5)' - character(len=*),parameter :: fmt2 = '(1x,a20," : ",f12.5)' - character(len=20) :: atmp - integer :: constraintype(8) + integer :: constraintype(8) write (iunit,'(1x,a)') '----------------' write (iunit,'(1x,a)') 'Calculation info' @@ -1066,6 +877,291 @@ subroutine calculation_info(self,iunit) return end subroutine calculation_info + + +!=========================================================================================! +!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< CALCULATION_SETTINGS associated routines +!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<-- check for missing settings in a calculation_settings object + subroutine calculation_settings_autocomplete(self,id) + implicit none + class(calculation_settings) :: self + integer,intent(in) :: id + integer :: i,j + character(len=50) :: nmbr + + !> add a short description + self%description = trim(jobdescription(self%id+1)) + call calculation_settings_shortflag(self) + + if (.not.allocated(self%calcspace)) then + !> I've decided to perform all calculations in a separate directory to + !> avoid accumulation of files in the main workspace + write (nmbr,'(i0)') id + self%calcspace = 'calculation.level.'//trim(nmbr) + end if + + if (self%pr .and. self%prch.ne.stdout) then + self%prch = self%prch+id + end if + end subroutine calculation_settings_autocomplete + +!>--- create a short calculation info flag + subroutine calculation_settings_shortflag(self) + implicit none + class(calculation_settings) :: self + integer :: i,j + + select case( self%id ) + case( jobtype%xtbsys ) + self%shortflag = 'xtb subprocess' + case( jobtype%generic ) + self%shortflag = 'generic subprocess' + case( jobtype%turbomole ) + self%shortflag = 'TURBOMOLE subprocess' + case( jobtype%orca ) + self%shortflag = 'ORCA subprocess' + case( jobtype%terachem ) + self%shortflag = 'TeraChem subprocess' + case( jobtype%tblite ) + select case (self%tblitelvl) + case (xtblvl%gfn2) + self%shortflag = 'GFN2-xTB' + case (xtblvl%gfn1) + self%shortflag = 'GFN1-xTB' + case (xtblvl%ipea1) + self%shortflag = 'IPEA1-xTB' + case (xtblvl%ceh) + self%shortflag = 'CEH' + case (xtblvl%eeq) + self%shortflag = 'EEQ(D4)' + case (xtblvl%param) + self%shortflag = 'parameter file: '//trim(self%tbliteparam) + end select + case( jobtype%gfn0 ) + self%shortflag = 'GFN0-xTB' + case( jobtype%gfn0occ ) + self%shortflag = 'GFN0-xTB*' + case( jobtype%gfnff ) + self%shortflag = 'GFN-FF' + case( jobtype%xhcff ) + self%shortflag = 'XHCFF' + case( jobtype%lj ) + self%shortflag = 'LJ' + case default + self%shortflag = 'undefined' + end select + if(allocated(self%solvmodel).and.allocated(self%solvent))then + self%shortflag = self%shortflag//'/'//trim(self%solvmodel) + self%shortflag = self%shortflag//'('//trim(self%solvent)//')' + endif + end subroutine calculation_settings_shortflag + +!>-- generate a unique print id for the calculation + subroutine calculation_settings_printid(self,thread,id) + implicit none + class(calculation_settings) :: self + integer,intent(in) :: thread,id + integer :: i,j,dum + character(len=50) :: nmbr + dum = 100*(thread+1) + dum = dum+id + self%prch = dum + end subroutine calculation_settings_printid + + subroutine calculation_dump_dipgrad(self,filename) + implicit none + class(calculation_settings) :: self + character(len=*),intent(in) :: filename + integer :: i,j,ich + if (.not.allocated(self%dipgrad)) return + open (newunit=ich,file=filename) + do i = 1,size(self%dipgrad,2) + write (ich,'(3f20.10)') self%dipgrad(1:3,i) + end do + close (ich) + end subroutine calculation_dump_dipgrad + +!=========================================================================================! + + subroutine calculation_settings_info(self,iunit) + implicit none + class(calculation_settings) :: self + integer,intent(in) :: iunit + integer :: i,j + character(len=*),parameter :: fmt1 = '(" :",2x,a20," : ",i0)' + character(len=*),parameter :: fmt2 = '(" :",2x,a20," : ",f0.5)' + character(len=*),parameter :: fmt3 = '(" :",2x,a20," : ",a)' + character(len=*),parameter :: fmt4 = '(" :",1x,a)' + character(len=20) :: atmp + + if (allocated(self%description)) then + write (iunit,'(" :",1x,a)') trim(self%description) + else + write (atmp,*) 'Job type' + write (iunit,fmt1) atmp,self%id + end if + !> more info + if (self%id == jobtype%tblite) then + select case (self%tblitelvl) + case (xtblvl%gfn2) + write (iunit,fmt4) 'GFN2-xTB level' + case (xtblvl%gfn1) + write (iunit,fmt4) 'GFN1-xTB level' + case (xtblvl%ceh) + write (iunit,fmt4) 'Charge Extended Hückel (CEH) model' + end select + end if + if (any((/jobtype%orca,jobtype%xtbsys,jobtype%turbomole, & + & jobtype%generic,jobtype%terachem/) == self%id)) then + write (iunit,'(" :",3x,a,a)') 'selected binary : ',trim(self%binary) + end if + if (self%refine_lvl > 0) then + write (atmp,*) 'refinement stage' + write (iunit,fmt1) atmp,self%refine_lvl + end if + + !> system data + write (atmp,*) 'Molecular charge' + write (iunit,fmt1) atmp,self%chrg + if (self%uhf /= 0) then + write (atmp,*) 'UHF parameter' + write (iunit,fmt1) atmp,self%uhf + end if + + if (allocated(self%solvmodel)) then + write (atmp,*) 'Solvation model' + write (iunit,fmt3) atmp,trim(self%solvmodel) + end if + if (allocated(self%solvent)) then + write (atmp,*) 'Solvent' + write (iunit,fmt3) atmp,trim(self%solvent) + end if + + !> xTB specific parameters + if (any((/jobtype%tblite,jobtype%xtbsys,jobtype%gfn0,jobtype%gfn0occ/) == self%id)) then + write (atmp,*) 'Fermi temperature' + write (iunit,fmt2) atmp,self%etemp + write (atmp,*) 'Accuracy' + write (iunit,fmt2) atmp,self%accuracy + write (atmp,*) 'max SCC cycles' + write (iunit,fmt1) atmp,self%maxscc + end if + + write (atmp,*) 'Reset data?' + if (self%apiclean) write (iunit,fmt3) atmp,'yes' + write (atmp,*) 'Read WBOs?' + if (self%rdwbo) write (iunit,fmt3) atmp,'yes' + write (atmp,*) 'Read dipoles?' + if (self%rddip) write (iunit,fmt3) atmp,'yes' + + if (self%ONIOM_highlowroot /= 0) then + select case (self%ONIOM_highlowroot) + case (1) + write (atmp,*) 'ONIOM frag ("high")' + case (2) + write (atmp,*) 'ONIOM frag ("low")' + case (3) + write (atmp,*) 'ONIOM frag ("root")' + end select + write (iunit,fmt1) trim(atmp),self%ONIOM_id + else + if(self%weight .ne. 1.0_wp)then + write (atmp,*) 'Weight' + write (iunit,fmt2) atmp,self%weight + endif + end if + + end subroutine calculation_settings_info + !=========================================================================================! subroutine create_calclevel_shortcut(self,levelstring) @@ -1101,68 +1197,10 @@ subroutine create_calclevel_shortcut(self,levelstring) self%id = jobtype%generic end select + call self%autocomplete(self%id) end subroutine create_calclevel_shortcut !=========================================================================================! - - subroutine calc_set_active(self,ids) - implicit none - class(calcdata) :: self - integer,intent(in) :: ids(:) - integer :: i,j,k,l - if (allocated(self%weightbackup)) deallocate (self%weightbackup) -!>--- on-the-fly multiscale definition - allocate (self%weightbackup(self%ncalculations),source=1.0_wp) - do i = 1,self%ncalculations -!>--- save backup weights - self%weightbackup(i) = self%calcs(i)%weight -!>--- set the weight of all unwanted calculations to 0 - if (.not.any(ids(:) .eq. i)) then - self%calcs(i)%weight = 0.0_wp - self%calcs(i)%active = .false. - else -!>--- and all other to 1 - self%calcs(i)%weight = 1.0_wp - self%calcs(i)%active = .true. - end if - end do - end subroutine calc_set_active - - subroutine calc_set_active_restore(self) - implicit none - class(calcdata) :: self - integer :: i,j,k,l - if (.not.allocated(self%weightbackup)) return - do i = 1,self%ncalculations -!>--- set all to active and restore saved weights - self%calcs(i)%weight = self%weightbackup(i) - self%calcs(i)%active = .true. - self%eweight(i) = self%weightbackup(i) - end do - deallocate (self%weightbackup) - end subroutine calc_set_active_restore - -!==========================================================================================! - - subroutine get_dipoles(calc,diptmp) -!********************************************* -!* Collect the x y and z dipole moments for -!* each calculation level in one array diptmp -!********************************************* - implicit none - type(calcdata),intent(inout) :: calc - real(wp),allocatable,intent(out) :: diptmp(:,:) - integer :: i,j,k,l,m - - m = calc%ncalculations - allocate (diptmp(3,m),source=0.0_wp) - do i = 1,m - if (calc%calcs(i)%rddip) then - diptmp(1:3,i) = calc%calcs(i)%dipole(1:3) - end if - end do - end subroutine get_dipoles - -!=========================================================================================! +!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ",a,i3,a:)') 'Calculation level ',k,'; '//calc%calcs(k)%shortflag - write (iunit,'(a12,a12,a10)') 'Atom A','Atom B','BO(A-B)' + write (iunit,'(a12,a12,a20)') 'Atom A','Atom B','BO(A-B)' do i = 1,mol%nat do j = 1,i-1 if (calc%calcs(k)%wbo(i,j) > 0.01_wp) then - write (iunit,*) i,j,calc%calcs(k)%wbo(i,j) + write (iunit,'(i12,i12,f20.8)') i,j,calc%calcs(k)%wbo(i,j) end if end do end do diff --git a/src/classes.f90 b/src/classes.f90 index e2dd4dde..7014d469 100644 --- a/src/classes.f90 +++ b/src/classes.f90 @@ -123,13 +123,14 @@ module crest_data !========================================================================================! !========================================================================================! - private - !========================================================================================! !========================================================================================! type :: constra +!**************************************************** +!* separate settings for LEGACY constraint handling +!**************************************************** integer :: ndim logical :: used = .false. logical :: NCI = .false. @@ -163,22 +164,29 @@ module crest_data !========================================================================================! type :: protobj - integer :: nfrag - integer :: newchrg - integer :: iter - real(wp) :: popthr - real(wp) :: ewin - integer :: swchrg !switch element charge - integer :: swat !switch element element - logical :: swelem !switch element to add to lmo lp pair? - logical :: allowFrag - logical :: threshsort !use ewin threshold - logical :: protdeprot !currently unused! - logical :: deprotprot !(tautomerize) do first deprotonation and then protonation - - logical :: strictPDT = .false. ! strict mode (i.e. bond constraints) for (de)protonation,tautomerization - logical :: fixPDT = .false. ! extension to the strict mode, fix heavy atom positions - logical :: ABcorrection = .false. +!************************************************************ +!* separate settings for protonation and related procedures +!************************************************************ + integer :: nfrag = 0 + integer :: newchrg = 0 + integer :: iter = 1 + real(wp) :: ewin = 30.0_wp !> separate EWIN threshold + integer :: swchrg = 1 !> switch element charge + integer :: swat = 1 !> switch element element + logical :: swelem = .false. !> switch element to add to lmo lp pair? + logical :: allowFrag = .false. !> allow fragmentation + logical :: threshsort = .false. !> use ewin threshold + logical :: protdeprot = .false. !> currently unused! + logical :: deprotprot = .false. !> (tautomerize) do first deprotonation and then protonation + + logical :: ffopt = .true. !> pre-optimize with a force-field to avoid high-energy artifacts + logical :: hnewopt = .true. !> optimization step only using the newly added atom + logical :: finalopt = .true. !> final optimization (only if ffopt is true) + + logical :: active_lmo(3) = (/.true.,.true.,.true./) !> consider pi, LP and delpi LMOs for protonation? + + logical :: strictPDT = .false. !> LEGACY: strict mode (i.e. bond constraints) for (de)protonation,tautomerization + logical :: fixPDT = .false. !> LEGACY: extension to the strict mode, fix heavy atom positions integer,allocatable :: atmap(:) @@ -190,6 +198,7 @@ module crest_data character(len=:),allocatable :: newligand !--- pka + logical :: ABcorrection = .false. integer :: h_acidic = 0 !which h atom to remove in pka script integer :: pka_mode = 0 !what to do in the pka calc. character(len=:),allocatable :: pka_baseinp !if a base file is read in instead @@ -205,8 +214,10 @@ module crest_data !========================================================================================! -!>--- ENTROPY mode setting object type :: entropyMTD +!********************************* +!* Separate ENTROPY mode settings +!********************************* integer :: nMDs ! number of static MTDs integer :: nBias ! number of Bias structures real(wp) :: nbiasgrow = 1.2d0 @@ -244,8 +255,11 @@ module crest_data end type entropyMTD !========================================================================================! -!>--- thermodynamics evaluation data + type :: thermodata +!***************************************** +!* separate thermodynamics evaluation data +!****************************************** real(wp) :: ithr = -50.0_wp !> imaginary mode inversion (in xtb -20.0) real(wp) :: fscal = 1.0_wp !> frequency scaling real(wp) :: sthr = 25.0_wp !> rot/vib interpol threshold (in xtb 50.0) @@ -265,8 +279,11 @@ module crest_data end type thermodata !========================================================================================! -!>--- storage of the reference (input structure) + type :: refdata +!****************************************** +!* separate storage of REFERENCE STRUCTURE +!****************************************** integer :: nat integer,allocatable :: at(:) real(wp),allocatable :: xyz(:,:) diff --git a/src/confparse.f90 b/src/confparse.f90 index f60f6bd2..19bd094a 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%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%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%protb%ffopt = .true. !> FF pre-optimization (CREST>3.0.2) env%pclean = .false. !> cleanup option for property mode !>--- options for principal component analysis (PCA) and clustering diff --git a/src/cregen.f90 b/src/cregen.f90 index 31d717b9..6004d5c5 100644 --- a/src/cregen.f90 +++ b/src/cregen.f90 @@ -23,12 +23,12 @@ !> This is a rewrite of the original routines since the old ones !> got a bit messy over time. !> the quickset variable can be used for some special runtypes: -!> quickset: 2 - do symmetry analysis -!> 3 - switch off equivalency analysis -!> 6 - energy sorting only -!> 9 - no sorting, only check groups -!> 12 - no topology check, turn ewin to infty -!> 13 - no topology check, ewin and rmsd checking (msreact settings) +!> quickset: 2 - do symmetry analysis +!> 3 - switch off equivalency analysis +!> 6,7 - energy sorting only with (7) or without (6) ewin energy cut-off +!> 9 - no sorting, only check groups +!> 12 - no topology check, turn ewin to infty +!> 13 - no topology check, ewin and rmsd checking (msreact settings) !=========================================================================================! !=========================================================================================! @@ -335,7 +335,7 @@ subroutine cregen_files(env,fname,oname,cname,simpleset,userinput,iounit) call remove(outfile) if (simpleset > 0) then select case (simpleset) - case (6,9,12) + case (6,7,9,12) iounit = stdout case default open (newunit=iounit,file=outfile) @@ -413,7 +413,7 @@ subroutine cregen_prout(env,simpleset,pr1,pr2,pr3,pr4) pr3 = .false. !> plain energy list pr4 = .false. !> group list printout - if (simpleset == 6) then + if (any(simpleset == (/6,7/))) then pr1 = .false. pr2 = .false. if (env%crestver .ne. crest_solv) pr3 = .true. @@ -489,7 +489,7 @@ subroutine cregen_director(env,simpleset,checkbroken,sortE,sortRMSD,sortRMSD2, & anal = .false. end if - if (simpleset == 6) then !energy sorting only + if (any(simpleset == (/6,7/))) then !energy sorting only checkbroken = .false. sortE = .true. sortRMSD = .false. @@ -575,7 +575,7 @@ subroutine cregen_filldata2(simpleset,ewin) implicit none integer,intent(in) :: simpleset real(wp),intent(out) :: ewin - if (simpleset == 6.or.simpleset == 12) then + if (any(simpleset == (/6,12/))) then ewin = huge(ewin) end if return @@ -1082,6 +1082,11 @@ subroutine cregen_esort(ch,nat,nall,xyz,comments,nallout,ewin) call stringqsort(nall,comments,1,nall,order) !>-- determine cut-off of energies + if (ewin < 9999.9_wp) then + write (ch,'('' sorting energy window (EWIN) :'',1x,f9.4,a)') ewin,' / kcal*mol⁻¹' + else + write (ch,'('' sorting energy window (EWIN) :'',3x,a,a)') '+∞',' / kcal*mol⁻¹' + end if emax = maxval(energies(:),1) de = (emax-energies(1))*autokcal if (de .gt. ewin) then @@ -1094,12 +1099,12 @@ subroutine cregen_esort(ch,nat,nall,xyz,comments,nallout,ewin) exit end if end do - write (ch,'('' number of removed by energy :'',i6)') (nall-nallout) - write (ch,'('' number of remaining points :'',i6)') nallout + write (ch,'('' number of removed by energy :'',3x,i0)') (nall-nallout) + write (ch,'('' number of remaining points :'',3x,i0)') nallout else nallout = nall end if - write (ch,*) 'reference state Etot :',energies(1) + write (ch,*) 'reference state Etot :',energies(1) deallocate (order,orderref) deallocate (energies) @@ -2670,12 +2675,6 @@ subroutine cregen_pr1(ch,env,nat,nall,rthr,bthr,pthr,ewin) write (ch,'('' RMSD threshold :'',f9.4)') rthr write (ch,'('' Bconst threshold :'',f9.4)') bthr write (ch,'('' population threshold :'',f9.4)') pthr - if (ewin < 9999.9_wp) then - write (ch,'('' conformer energy window /kcal :'',f9.4)') ewin - else - write (ch,'('' conformer energy window /kcal :'',6x,a)') '+∞' - end if - return end subroutine cregen_pr1 @@ -2924,15 +2923,16 @@ subroutine cregen_pr3(ch,infile,nall,comments) end do write (ch,*) - write (ch,'(a)') '===================================================' - write (ch,'(a)') '============= ordered structure list ==============' - write (ch,'(a)') '===================================================' - write (ch,'(a,a,a)') ' written to file <',trim(infile),'>' + write (ch,'(a)') '=====================================================' + write (ch,'(a)') '============== ordered structure list ===============' + write (ch,'(a)') '=====================================================' + write (ch,'(a,a,a)') ' written to file <',trim(infile),'>' write (ch,*) - write (ch,'('' structure ΔE(kcal/mol) Etot(Eh)'')') + write (ch,'(a10,4x,a15,a25)') 'structure','ΔE(kcal/mol)','Etot(Eh)' + !write (ch,'('' structure ΔE(kcal/mol) Etot(Eh)'')') do i = 1,nall dE = (er(i)-er(1))*autokcal - write (ch,'(i10,3x,F12.2,2x,F14.6)') i,dE,er(i) + write (ch,'(i10,3x,F15.4,F25.10)') i,dE,er(i) end do write (ch,*) diff --git a/src/legacy_wrappers.f90 b/src/legacy_wrappers.f90 index b0ddd427..daf9a909 100644 --- a/src/legacy_wrappers.f90 +++ b/src/legacy_wrappers.f90 @@ -50,7 +50,7 @@ subroutine env2calc(env,calc,molin) cal%chrg = env%chrg !>-- obtain WBOs OFF by default cal%rdwbo = .false. - cal%rddip = .true. + cal%rddip = .false. !> except for SP runtype (from command line!) if (env%crestver == crest_sp) then cal%rdwbo = .true. @@ -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 !========================================================================================! diff --git a/src/minitools.f90 b/src/minitools.f90 index eb41198d..a6be69f2 100644 --- a/src/minitools.f90 +++ b/src/minitools.f90 @@ -25,7 +25,7 @@ subroutine splitfile(fname,up,low) !******************************************************** use crest_parameters use iomod - use strucrd,only:rdensembleparam,rdensemble,wrxyz + use strucrd,only: rdensemble,coord implicit none character(len=*) :: fname integer :: up,low @@ -38,6 +38,7 @@ subroutine splitfile(fname,up,low) integer,allocatable :: at(:) integer :: i,r logical :: ex + type(coord),allocatable :: structures(:) inquire (file=fname,exist=ex) if (.not.ex) then @@ -47,9 +48,10 @@ subroutine splitfile(fname,up,low) call getcwd(thispath) !current dir= thispath - call rdensembleparam(fname,nat,nall) - allocate (xyz(3,nat,nall),at(nat)) - call rdensemble(fname,nat,nall,at,xyz) + !call rdensembleparam(fname,nat,nall) + !allocate (xyz(3,nat,nall),at(nat)) + !call rdensemble(fname,nat,nall,at,xyz) + call rdensemble(fname,nall,structures) r = makedir("SPLIT") !create new directory call chdir("SPLIT") @@ -73,7 +75,8 @@ subroutine splitfile(fname,up,low) write (tmppath2,'(a,i0)') "STRUC",i r = makedir(trim(tmppath2)) call chdir(tmppath2) - call wrxyz("struc.xyz",nat,at,xyz(:,:,i)) + !call wrxyz("struc.xyz",nat,at,xyz(:,:,i)) + call structures(i)%write("struc.xyz") call chdir(tmppath1) end do diff --git a/src/optimize/optimize_module.f90 b/src/optimize/optimize_module.f90 index 431da6b7..42fbe67f 100644 --- a/src/optimize/optimize_module.f90 +++ b/src/optimize/optimize_module.f90 @@ -65,6 +65,13 @@ subroutine optimize_geometry(mol,molnew,calc,etot,grd,pr,wr,iostatus) molnew%nat = mol%nat !$omp end critical + !> Check for optimization-individual calculation setup + if(calc%optnewinit)then + !$omp critical + call calc%dealloc_params() + !$omp end critical + endif + !> initial singlepoint call engrad(molnew,calc,etot,grd,iostatus) @@ -85,7 +92,7 @@ subroutine optimize_geometry(mol,molnew,calc,etot,grd,pr,wr,iostatus) write(stdout,'(a)') 'Unknown optimization engine!' stop end select - + molnew%energy = etot return end subroutine optimize_geometry @@ -132,6 +139,10 @@ subroutine print_opt_data(calc,ich) write (ich,'(1x,a,e10.3,a,e10.3,a)') 'E/G convergence criteria: ',& & ethr,' Eh,',gthr,' Eh/a0' + write (ich,'(1x,a,i0)') 'maximum optimization steps: ',calc%maxcycle + end subroutine print_opt_data + +!========================================================================================! !========================================================================================! end module optimize_module diff --git a/src/parsing/confparse2.f90 b/src/parsing/confparse2.f90 index 108bff39..2dcc46f3 100644 --- a/src/parsing/confparse2.f90 +++ b/src/parsing/confparse2.f90 @@ -62,6 +62,7 @@ subroutine parseinputfile(env,fname) type(mddata) :: mddat logical :: ex,l1,l2 integer :: i,j,k,l + integer :: readstatus !>--- check for the input files' existence inquire (file=fname,exist=ex) @@ -74,12 +75,15 @@ subroutine parseinputfile(env,fname) !>--- read the file into the object 'dict' call parse_input(fname,dict) - call dict%print() + call dict%print2() + +!>--- sanity check for input files + readstatus = 0 !> has to remain 0, or something went wrong !>--- parse all root-level key-value pairs do i = 1,dict%nkv kv = dict%kv_list(i) - call parse_main_auto(env,kv) + call parse_main_auto(env,kv,readstatus) end do !>------------------------------------------------------ !> After this point I assume an input structure was @@ -89,12 +93,12 @@ subroutine parseinputfile(env,fname) !>--- parse all objects that write to env or global data do i = 1,dict%nblk blk = dict%blk_list(i) - call parse_main_blk(env,blk) + call parse_main_blk(env,blk,readstatus) end do !>--- check objects for a calculation setup ! i.e., all [calculation] and [[calculation.*]] blocks - call parse_calculation_data(env,newcalc,dict,l1) + call parse_calculation_data(env,newcalc,dict,l1,readstatus) if (l1) then env%calc = newcalc call env_calcdat_specialcases(env) @@ -102,11 +106,17 @@ subroutine parseinputfile(env,fname) !>--- check for molecular dynamics setup ! i.e., all [dynamics] and [[dynamics.*]] blocks - call parse_dynamics_data(env,mddat,dict,l1) + call parse_dynamics_data(env,mddat,dict,l1,readstatus) if (l1) then env%mddat = mddat end if +!>--- terminate if there were any unrecognized keywords + if(readstatus /= 0)then + write(stdout, '(i0,a)') readstatus,' error(s) while reading input file' + error stop + endif + !>--- check for lwONIOM setup (will be read at end of confparse) do i = 1,dict%nblk if (dict%blk_list(i)%header == 'lwoniom') then diff --git a/src/parsing/parse_block.f90 b/src/parsing/parse_block.f90 index 110a876a..44678796 100644 --- a/src/parsing/parse_block.f90 +++ b/src/parsing/parse_block.f90 @@ -15,6 +15,7 @@ module parse_block contains procedure :: addkv => blk_addkv procedure :: print => blk_print + procedure :: print2 => blk_print2 procedure :: fmt_header => blk_fmt_header procedure :: deallocate => blk_deallocate end type datablock @@ -95,6 +96,20 @@ subroutine blk_print(self) end if end do end subroutine blk_print + + subroutine blk_print2(self) + class(datablock) :: self + integer :: i + if(index(self%header,'.').ne.0)then + write (stdout,'("*",1x,a)') '[['//self%header//']]' + else + write (stdout,'("*",1x,a)') '['//self%header//']' + endif + do i = 1,self%nkv + write (stdout,'("*",1x,a)') trim(self%kv_list(i)%print2()) + end do + end subroutine blk_print2 + !========================================================================================! !> the following routines are only used in the fallback parsing routines @@ -197,7 +212,7 @@ subroutine deallocate_block(self) end subroutine deallocate_block !========================================================================================! -!> deallocate block data +!> print block data subroutine print_block(self) implicit none class(parseblock) :: self @@ -215,6 +230,25 @@ subroutine print_block(self) return end subroutine print_block + subroutine print_block2(self) + implicit none + class(parseblock) :: self + integer :: i + + write (*,*) + if (allocated(self%header)) then + write (*,*) trim(self%header) + if (allocated(self%content)) then + do i = 1,self%len + write (*,*) self%content(i) + end do + end if + end if + return + end subroutine print_block2 + + + !=======================================================================================! !> check if string is a toml block header function isheader(str) diff --git a/src/parsing/parse_calcdata.f90 b/src/parsing/parse_calcdata.f90 index 977ae89b..0f93444b 100644 --- a/src/parsing/parse_calcdata.f90 +++ b/src/parsing/parse_calcdata.f90 @@ -43,49 +43,36 @@ module parse_calcdata !>-- routines for parsing a calculation_settings object interface parse_setting module procedure :: parse_setting_auto - module procedure :: parse_setting_float - module procedure :: parse_setting_int - module procedure :: parse_setting_c - module procedure :: parse_setting_bool end interface parse_setting !>-- routines for parsing a calcdata object interface parse_calc module procedure :: parse_calc_auto - module procedure :: parse_calc_float - module procedure :: parse_calc_int - module procedure :: parse_calc_c - module procedure :: parse_calc_bool end interface parse_calc !>-- routines for parsing a mddata object interface parse_md module procedure :: parse_md_auto - module procedure :: parse_md_float - module procedure :: parse_md_int - module procedure :: parse_md_c - module procedure :: parse_md_bool end interface parse_md !>-- routines for parsing a mtdpot object interface parse_mtd module procedure :: parse_metadyn_auto - module procedure :: parse_mtd_float - module procedure :: parse_mtd_int - module procedure :: parse_mtd_c - module procedure :: parse_mtd_bool end interface parse_mtd public :: parse_calculation_data public :: parse_dynamics_data + character(len=*),parameter,private :: fmturk = '("unrecognized KEYWORD in ",a," : ",a)' + character(len=*),parameter,private :: fmtura = '("unrecognized ARGUMENT : ",a)' + !========================================================================================! !========================================================================================! contains !> MODULE PROCEDURES START HERE !========================================================================================! !========================================================================================! - subroutine parse_calculation_data(env,calc,dict,included) + subroutine parse_calculation_data(env,calc,dict,included,istat) implicit none type(systemdata) :: env type(calcdata) :: calc @@ -94,23 +81,24 @@ subroutine parse_calculation_data(env,calc,dict,included) type(calculation_settings) :: newjob,newjob2 type(constraint) :: newcstr integer :: i,j,k,l - logical,intent(out) :: included + logical,intent(out) :: included + integer,intent(inout) :: istat type(coord) :: moltmp included = .false. call calc%reset() call env%ref%to(moltmp) - call axis(moltmp%nat,moltmp%at,moltmp%xyz) + call axis(moltmp%nat,moltmp%at,moltmp%xyz) do i = 1,dict%nblk call blk%deallocate() blk = dict%blk_list(i) if (blk%header == 'calculation') then included = .true. - call parse_calcdat(env,blk,calc) + call parse_calcdat(env,blk,calc,istat) else if (blk%header == 'calculation.level') then - call parse_leveldata(env,blk,newjob) + call parse_leveldata(env,blk,newjob,istat) call newjob%autocomplete(calc%ncalculations+1) call calc%add(newjob) included = .true. @@ -120,9 +108,9 @@ subroutine parse_calculation_data(env,calc,dict,included) if (allocated(calc%calcs)) deallocate (calc%calcs) calc%ncalculations = 0 calc%id = -1 - call parse_leveldata(env,blk,newjob) + call parse_leveldata(env,blk,newjob,istat) !>-- S0 setup - call parse_leveldata(env,blk,newjob) + call parse_leveldata(env,blk,newjob,istat) newjob%uhf = 0 newjob%calcspace = 's0' call calc%add(newjob) @@ -132,12 +120,12 @@ subroutine parse_calculation_data(env,calc,dict,included) call calc%add(newjob) included = .true. - else if (index(blk%header,'calculation.constraint').ne.0) then - call parse_constraintdat(env,moltmp,blk,calc) + else if (index(blk%header,'calculation.constraint') .ne. 0) then + call parse_constraintdat(env,moltmp,blk,calc,istat) included = .true. else if (blk%header == 'calculation.scans') then - call parse_scandat(blk,calc) + call parse_scandat(blk,calc,istat) included = .true. end if @@ -150,7 +138,7 @@ end subroutine parse_calculation_data !========================================================================================! - subroutine parse_leveldata(env,blk,job) + subroutine parse_leveldata(env,blk,job,istat) !********************************************************** !* The following routines are used to !* read information into the "calculation_settings" object @@ -159,6 +147,8 @@ subroutine parse_leveldata(env,blk,job) type(systemdata),intent(inout) :: env type(datablock),intent(in) :: blk type(calculation_settings),intent(out) :: job + integer,intent(inout) :: istat + logical :: rd integer :: i call job%deallocate() if ((blk%header .ne. 'calculation.level').and. & @@ -166,91 +156,55 @@ subroutine parse_leveldata(env,blk,job) return end if do i = 1,blk%nkv - call parse_setting(env,job,blk%kv_list(i)) + call parse_setting_auto(env,job,blk%kv_list(i),rd) + if (.not.rd) then + istat = istat+1 + write (stdout,fmturk) '[['//blk%header//']]-block',blk%kv_list(i)%key + end if end do return end subroutine parse_leveldata - subroutine parse_setting_auto(env,job,kv) + subroutine parse_setting_auto(env,job,kv,rd) implicit none type(systemdata),intent(inout) :: env type(calculation_settings) :: job type(keyvalue) :: kv - !> first, go through settings with fixed type expactations - select case (kv%id) - case (1) !> float - call parse_setting(job,kv%key,kv%value_f) - case (2) !> int - call parse_setting(job,kv%key,kv%value_i) - case (3) !> bool - call parse_setting(job,kv%key,kv%value_b) - case (4) !> string - call parse_setting(env,job,kv%key,kv%value_c) - case (6,7) !> int/float array - call parse_setting_array(job,kv,kv%key) - end select - !> then, all others by key, automatic + logical,intent(out) :: rd + logical :: ex + rd = .true. select case (kv%key) - case default - continue - end select - end subroutine parse_setting_auto - subroutine parse_setting_float(job,key,val) - implicit none - type(calculation_settings) :: job - character(len=*) :: key - real(wp) :: val - select case (key) - case ('uhf') - job%uhf = nint(val) - case ('chrg','charge') - job%chrg = nint(val) + +!>--- floats case ('etemp') - job%etemp = val + job%etemp = kv%value_f case ('accuracy') - job%accuracy = val + job%accuracy = kv%value_f case ('weight') - job%weight = val + job%weight = kv%value_f case ('pressure') - job%extpressure = val + job%extpressure = kv%value_f case ('proberad') - job%proberad = val - end select - return - end subroutine parse_setting_float - subroutine parse_setting_int(job,key,val) - implicit none - type(calculation_settings) :: job - character(len=*) :: key - integer :: val - select case (key) + job%proberad = kv%value_f + +!>--- integers case ('uhf','multiplicity') - job%uhf = val + job%uhf = kv%value_i case ('chrg','charge') - job%chrg = val + job%chrg = kv%value_i case ('id') - job%id = val + job%id = kv%value_i case ('maxscc') - job%maxscc = val - case ('tblite_level','tblite_hamiltonian') - job%tblitelvl = val + job%maxscc = kv%value_i case ('lebedev') - job%ngrid = val + job%ngrid = kv%value_i case ('vdwset') - job%vdwset = val - end select - return - end subroutine parse_setting_int - subroutine parse_setting_c(env,job,key,val) - implicit none - type(systemdata),intent(inout) :: env - type(calculation_settings) :: job - character(len=*) :: key - character(len=*) :: val - logical :: ex - select case (key) + job%vdwset = kv%value_i + case ('config') + call job%addconfig(kv%value_ia) +!>--- strings case ('method') - select case (val) + select case (kv%value_c) case ('gfn-xtb','gfn','xtb') job%id = jobtype%xtbsys case ('generic') @@ -289,26 +243,29 @@ subroutine parse_setting_c(env,job,key,val) job%id = jobtype%lj case default job%id = jobtype%unknown + !>--- keyword was recognized, but invalid argument supplied + write (stdout,fmtura) kv%value_c + error stop end select case ('bin','binary','script') - job%binary = val + job%binary = kv%value_c case ('flags') - job%other = val + job%other = kv%value_c !> don't. !case ('sys','syscall','systemcall') ! job%systemcall = val case ('calcspace','dir') - job%calcspace = val + job%calcspace = kv%value_c case ('gradfile') - job%gradfile = val + job%gradfile = kv%value_c case ('gradtype') - select case (val) + select case (kv%value_c) case ('engrad','xtb','orca') job%gradtype = gradtype%engrad case ('turbomole','tm') @@ -317,19 +274,22 @@ subroutine parse_setting_c(env,job,key,val) job%gradtype = gradtype%unknown case default job%gradtype = gradtype%unknown + !>--- keyword was recognized, but invalid argument supplied + write (stdout,fmtura) kv%value_c + error stop end select case ('gradkey') - job%gradkey = val + job%gradkey = kv%value_c case ('gradmt') - job%gradfmt = conv2gradfmt(val) + job%gradfmt = conv2gradfmt(kv%value_c) case ('efile') - job%efile = val + job%efile = kv%value_c case ('tblite_level','tblite_hamiltonian') - select case (val) + select case (kv%value_c) case ('gfn2','gfn2-xtb') job%tblitelvl = xtblvl%gfn2 case ('gfn1','gfn1-xtb') @@ -341,28 +301,31 @@ subroutine parse_setting_c(env,job,key,val) job%rdgrad = .false. case ('eeq','d4eeq') job%tblitelvl = xtblvl%eeq - job%rdgrad=.false. + job%rdgrad = .false. case default job%tblitelvl = xtblvl%unknown + !>--- keyword was recognized, but invalid argument supplied + write (stdout,fmtura) kv%value_c + error stop end select - case('tblite_param') - job%tbliteparam = val + case ('tblite_param') + job%tbliteparam = kv%value_c job%tblitelvl = xtblvl%param case ('orca_cmd') job%id = jobtype%orca - job%ORCA%cmd = val - job%binary = val + job%ORCA%cmd = kv%value_c + job%binary = kv%value_c case ('orca_template') job%id = jobtype%orca - call job%ORCA%read(val) + call job%ORCA%read(kv%value_c) case ('gbsa','alpb','cpcm') - job%solvmodel = key - job%solvent = val + job%solvmodel = kv%key + job%solvent = kv%value_c case ('refine','refinement') - select case (val) + select case (kv%value_c) case ('sp','singlepoint') job%refine_lvl = refine%singlepoint case ('add','correction') @@ -371,105 +334,94 @@ subroutine parse_setting_c(env,job,key,val) job%refine_lvl = refine%geoopt case default job%refine_lvl = refine%non + !>--- keyword was recognized, but invalid argument supplied + write (stdout,fmtura) kv%value_c + error stop end select - case('restartfile','topo','reftopo') - inquire(file=val,exist=ex) - if(ex)then + case ('restartfile','topo','reftopo') + inquire (file=kv%value_c,exist=ex) + if (ex) then job%restart = .true. - job%restartfile = val + job%restartfile = kv%value_c else - write(stderr,'(a,a,a)') 'specified restart file ',val,' does not exist' + write (stderr,'(a,a,a)') 'specified restart file ',kv%value_c,' does not exist' error stop - endif - case('refgeo','refxyz') - inquire(file=val,exist=ex) - if(ex)then - job%refgeo = val + end if + case ('refgeo','refxyz') + inquire (file=kv%value_c,exist=ex) + if (ex) then + job%refgeo = kv%value_c else - write(stderr,'(a,a,a)') 'specified reference geometry file ',val,' does not exist' + write (stderr,'(a,a,a)') 'specified reference geometry file ',kv%value_c,' does not exist' error stop - endif - case('parametrisation') - inquire(file=val,exist=ex) - if(ex)then - job%parametrisation = val + end if + case ('parametrisation') + inquire (file=kv%value_c,exist=ex) + if (ex) then + job%parametrisation = kv%value_c else - write(stderr,'(a,a,a)') 'specified parametrisation file ',val,' does not exist' + write (stderr,'(a,a,a)') 'specified parametrisation file ',kv%value_c,' does not exist' error stop - endif - case('refchrg','refcharges') - inquire(file=val,exist=ex) - if(ex)then - job%refcharges = val + end if + case ('refchrg','refcharges') + inquire (file=kv%value_c,exist=ex) + if (ex) then + job%refcharges = kv%value_c else - write(stderr,'(a,a,a)') 'specified reference charge file ',val,' does not exist' + write (stderr,'(a,a,a)') 'specified reference charge file ',kv%value_c,' does not exist' error stop - endif - + end if case ('print') - select case (val) - case ('true','yes') - job%pr = .true. - case ('false','no') - job%pr = .false. - case ('append','cont','continuous') - job%pr = .true. - job%prappend = .true. + select case (kv%id) + case (2) + select case (kv%value_c) + case ('true','yes') + job%pr = .true. + case ('false','no') + job%pr = .false. + case ('append','cont','continuous') + job%pr = .true. + job%prappend = .true. + end select + case (3) + job%pr = kv%value_b end select if (job%pr) job%prch = 999 !> the actual ID will be generated automatically - case('getsasa') - call get_atlist(env%ref%nat,job%getsasa,val,env%ref%at) + case ('getsasa') + call get_atlist(env%ref%nat,job%getsasa,kv%value_c,env%ref%at) - end select - return - end subroutine parse_setting_c - subroutine parse_setting_bool(job,key,val) - implicit none - type(calculation_settings) :: job - character(len=*) :: key - logical :: val - select case (key) +!>--- booleans case ('rdwbo') - job%rdwbo = val + job%rdwbo = kv%value_b case ('rddip','rddipole') - job%rddip = val + job%rddip = kv%value_b case ('rdqat','rdchrg') - job%rdqat = val + job%rdqat = kv%value_b case ('dumpq','dumpchrg') - job%rdqat = val - job%dumpq = val + job%rdqat = kv%value_b + job%dumpq = kv%value_b case ('dipgrad') - job%rddipgrad = val + job%rddipgrad = kv%value_b case ('rdgrad') - job%rdgrad = val + job%rdgrad = kv%value_b case ('refresh') - job%apiclean = val + job%apiclean = kv%value_b case ('lmo','lmocent') - job%getlmocent = val - case ('print') - job%pr = val - if (val) job%prch = 999 !> the actual ID will be generated automatically - end select - return - end subroutine parse_setting_bool - subroutine parse_setting_array(job,kv,key) - implicit none - type(calculation_settings) :: job - type(keyvalue) :: kv - character(len=*) :: key - select case (key) - case ('config') - call job%addconfig(kv%value_ia) + job%getlmocent = kv%value_b + + case default + !>--- keyword not correctly read/found + rd = .false. + continue end select - return - end subroutine parse_setting_array + end subroutine parse_setting_auto !========================================================================================! - subroutine parse_calcdat(env,blk,calc) + subroutine parse_calcdat(env,blk,calc,istat) !*********************************************** !* The following routines are used to !* read information into the "calcdata" object @@ -478,100 +430,72 @@ subroutine parse_calcdat(env,blk,calc) type(systemdata),intent(inout) :: env type(datablock),intent(in) :: blk type(calcdata),intent(inout) :: calc + integer,intent(inout) :: istat integer :: i + logical :: rd if (blk%header .ne. 'calculation') return do i = 1,blk%nkv - call parse_calc(env,calc,blk%kv_list(i)) + call parse_calc_auto(env,calc,blk%kv_list(i),rd) + if (.not.rd) then + istat = istat+1 + write (stdout,fmturk) '[calculation]-block',blk%kv_list(i)%key + end if end do return end subroutine parse_calcdat - subroutine parse_calc_auto(env,calc,kv) + subroutine parse_calc_auto(env,calc,kv,rd) implicit none type(systemdata),intent(inout) :: env type(calcdata) :: calc type(keyvalue) :: kv - select case (kv%id) - case (1) !> float - call parse_calc(env,calc,kv%key,kv%value_f) - case (2) !> int - call parse_calc(env,calc,kv%key,kv%value_i) - case (3) !> bool - call parse_calc(env,calc,kv%key,kv%value_b) - case (4) !> string - call parse_calc(env,calc,kv%key,kv%value_c) - end select - !> other, with multiple or raw type + logical,intent(out) :: rd + logical,allocatable :: atlist(:) + rd = .true. select case (kv%key) case ('optlev','ancopt_level') env%optlev = optlevnum(kv%rawvalue) - end select - end subroutine parse_calc_auto - subroutine parse_calc_float(env,calc,key,val) - implicit none - type(systemdata),intent(inout) :: env - type(calcdata) :: calc - character(len=*) :: key - real(wp) :: val - select case (key) + +!>--- floats case ('converge_e','ethr_opt') - calc%ethr_opt = val !> optimization ΔE convergenve threshold (Ha) + calc%ethr_opt = kv%value_f !> optimization ΔE convergenve threshold (Ha) case ('converge_g','gthr_opt','rmsforce') - calc%gthr_opt = val !> optimization RMS convergence threshold (Ha/a0) + calc%gthr_opt = kv%value_f !> optimization RMS convergence threshold (Ha/a0) case ('maxerise') - calc%maxerise = val !> optimization max E rise (Ha) + calc%maxerise = kv%value_f !> optimization max E rise (Ha) case ('displ_opt','maxdispl') - calc%maxdispl_opt = val !> optimization step size/scaling + calc%maxdispl_opt = kv%value_f !> optimization step size/scaling case ('hguess') - calc%hguess = val !> guess for the initial hessian - - case default - return - end select - return - end subroutine parse_calc_float - subroutine parse_calc_int(env,calc,key,val) - implicit none - type(systemdata),intent(inout) :: env - type(calcdata) :: calc - character(len=*) :: key - integer :: val - select case (key) - case ('id','type') - calc%id = val + calc%hguess = kv%value_f !> guess for the initial hessian +!>--- integers case ('maxcycle') - calc%maxcycle = val !> optimization max cycles + calc%maxcycle = kv%value_i !> optimization max cycles - case default - return - end select - return - end subroutine parse_calc_int - subroutine parse_calc_c(env,calc,key,val) - implicit none - type(systemdata),intent(inout) :: env - type(calcdata) :: calc - character(len=*) :: key - character(len=*) :: val - logical,allocatable :: atlist(:) - select case (key) - case ('type') - select case (val) - case ('mecp') - calc%id = -1 - case default - calc%id = 0 +!>--- strings + case ('id','type') + !> (OLD setting) calculation type + select case (kv%id) + case (2) + calc%id = kv%value_i + case (4) + select case (kv%value_c) + case ('mecp') + calc%id = -1 + case default + calc%id = 0 + end select end select + case ('elog') - calc%elog = val + calc%elog = kv%value_c calc%pr_energies = .true. case ('hess_update','hupdate') - select case (val) !> Hessian updates in geom. Opt. + select case (kv%value_c) !> Hessian updates in geom. Opt. case ('bfgs') calc%iupdat = 0 case ('powell') @@ -583,11 +507,13 @@ subroutine parse_calc_c(env,calc,key,val) case ('schlegel') calc%iupdat = 4 case default - calc%iupdat = 0 + !>--- keyword was recognized, but invalid argument supplied + write (stdout,fmtura) kv%value_c + error stop end select case ('opt','opt_engine','opt_algo') - select case (val) + select case (kv%value_c) case ('ancopt','rfo-anc') calc%opt_engine = 0 case ('lbfgs','l-bfgs') @@ -596,40 +522,32 @@ subroutine parse_calc_c(env,calc,key,val) calc%opt_engine = 2 case ('gd','gradient descent') calc%opt_engine = -1 + case default + !>--- keyword was recognized, but invalid argument supplied + write (stdout,fmtura) kv%value_c + error stop end select case ('freeze') - call get_atlist(env%ref%nat,atlist,val,env%ref%at) + call get_atlist(env%ref%nat,atlist,kv%value_c,env%ref%at) calc%nfreeze = count(atlist) call move_alloc(atlist,calc%freezelist) - case default - return - end select - return - end subroutine parse_calc_c - subroutine parse_calc_bool(env,calc,key,val) - implicit none - type(systemdata),intent(inout) :: env - type(calcdata) :: calc - character(len=*) :: key - logical :: val - select case (key) +!>--- booleans case ('eprint') - calc%pr_energies = val + calc%pr_energies = kv%value_b case ('exact_rf') - calc%exact_rf = val + calc%exact_rf = kv%value_b case default - return + rd = .false. end select - return - end subroutine parse_calc_bool + end subroutine parse_calc_auto !========================================================================================! - subroutine parse_constraintdat(env,mol,blk,calc) + subroutine parse_constraintdat(env,mol,blk,calc,istat) !************************************************* !* The following routines are used to !* read information into the "constraint" object @@ -640,24 +558,21 @@ subroutine parse_constraintdat(env,mol,blk,calc) type(coord),intent(inout) :: mol type(datablock),intent(in) :: blk type(calcdata),intent(inout) :: calc + integer,intent(inout) :: istat logical :: success type(constraint) :: constr integer :: i - !type(coord) :: mol - logical,allocatable :: atlist(:) - if (blk%header .ne. 'calculation.constraints' .and. & + logical :: rd + if (blk%header .ne. 'calculation.constraints'.and. & & blk%header .ne. 'calculation.constraint') return success = .false. call constr%deallocate() do i = 1,blk%nkv - call parse_constraint_auto(env,constr,blk%kv_list(i),success) - - select case (blk%kv_list(i)%key) - case ('freeze') - call get_atlist(env%ref%nat,atlist,blk%kv_list(i)%rawvalue,env%ref%at) - calc%nfreeze = count(atlist) - call move_alloc(atlist,calc%freezelist) - end select + call parse_constraint_auto(env,calc,constr,blk%kv_list(i),success,rd) + if (.not.rd) then + istat = istat+1 + write (stdout,fmturk) '[['//blk%header//']]-block',blk%kv_list(i)%key + end if end do if (success) then call constr%complete(mol) @@ -665,18 +580,26 @@ subroutine parse_constraintdat(env,mol,blk,calc) end if return end subroutine parse_constraintdat - subroutine parse_constraint_auto(env,constr,kv,success) + subroutine parse_constraint_auto(env,calc,constr,kv,success,rd) implicit none type(systemdata) :: env type(keyvalue) :: kv type(constraint) :: constr logical,intent(inout) :: success + type(calcdata),intent(inout) :: calc real(wp) :: dum1,dum2,dum3,dum4 real(wp) :: rabc(3) integer :: atm1,atm2,atm3,atm4,n,k,j logical,allocatable :: atlist(:) - !success = .false. + logical,intent(out) :: rd + rd = .true. select case (kv%key) + + case ('freeze') + call get_atlist(env%ref%nat,atlist,kv%rawvalue,env%ref%at) + calc%nfreeze = count(atlist) + call move_alloc(atlist,calc%freezelist) + case ('type') !> the type of constraint select case (kv%value_c) case ('bond','bonds'); constr%type = 1 @@ -688,6 +611,10 @@ subroutine parse_constraint_auto(env,constr,kv,success) case ('bondrange'); constr%type = 8 case ('gapdiff'); constr%type = -1 case ('gapdiff2','mecp'); constr%type = -2 + case default + !>--- keyword was recognized, but invalid argument supplied + write (stdout,fmtura) kv%value_c + error stop end select if (constr%type /= 0) success = .true. @@ -728,7 +655,7 @@ subroutine parse_constraint_auto(env,constr,kv,success) constr%atms(k) = j end if end do - deallocate (atlist) + deallocate (atlist) end if constr%n = n @@ -752,11 +679,11 @@ subroutine parse_constraint_auto(env,constr,kv,success) end select case ('wscal') !> scaling factor if the wall potential is automatically set up - if(kv%id == valuetypes%int)then - constr%wscal = max(0.0_wp, real(kv%value_i)) - elseif(kv%id == valuetypes%float)then - constr%wscal = max(0.0_wp, kv%value_f) - endif + if (kv%id == valuetypes%int) then + constr%wscal = max(0.0_wp,real(kv%value_i)) + elseif (kv%id == valuetypes%float) then + constr%wscal = max(0.0_wp,kv%value_f) + end if !>--- the following are for specifiying keywords in a single line !>--- I don't know it was wise to code them like this because it's hacky, @@ -873,6 +800,7 @@ subroutine parse_constraint_auto(env,constr,kv,success) !>-------------- !>-------------- case default + rd = .false. return end select @@ -881,7 +809,7 @@ end subroutine parse_constraint_auto !========================================================================================! - subroutine parse_scandat(blk,calc) + subroutine parse_scandat(blk,calc,istat) !******************************************* !* The following routines are used to !* read information into the "scan" object @@ -890,27 +818,34 @@ subroutine parse_scandat(blk,calc) implicit none type(datablock),intent(in) :: blk type(calcdata),intent(inout) :: calc + integer,intent(inout) :: istat logical :: success type(scantype) :: scn integer :: i + logical :: rd if (blk%header .ne. 'calculation.scans') return do i = 1,blk%nkv - call parse_scan_auto(scn,blk%kv_list(i),success) + call parse_scan_auto(scn,blk%kv_list(i),success,rd) + if (.not.rd) then + istat = istat+1 + write (stdout,fmturk) '[['//blk%header//']]-block',blk%kv_list(i)%key + end if if (success) then call calc%add(scn) end if end do return end subroutine parse_scandat - subroutine parse_scan_auto(scn,kv,success) + subroutine parse_scan_auto(scn,kv,success,rd) implicit none type(keyvalue) :: kv type(scantype) :: scn - logical,intent(out) :: success + logical,intent(out) :: success,rd real(wp) :: dum1,dum2,dum3 integer :: atm1,atm2,atm3,atm4 integer :: nsteps success = .false. + rd = .true. call scn%deallocate() select case (kv%key) case ('bond','distance') @@ -974,6 +909,7 @@ subroutine parse_scan_auto(scn,kv,success) end if case default + rd = .false. return end select @@ -982,7 +918,7 @@ end subroutine parse_scan_auto !========================================================================================! - subroutine parse_dynamics_data(env,mddat,dict,included) + subroutine parse_dynamics_data(env,mddat,dict,included,istat) !********************************************* !* The following routines are used to !* read information into the "mddata" object @@ -996,6 +932,7 @@ subroutine parse_dynamics_data(env,mddat,dict,included) type(constraint) :: newcstr integer :: i,j,k,l logical,intent(out) :: included + integer,intent(inout) :: istat included = .false. @@ -1004,9 +941,9 @@ subroutine parse_dynamics_data(env,mddat,dict,included) blk = dict%blk_list(i) if (blk%header == 'dynamics') then included = .true. - call parse_mddat(env,blk,mddat) + call parse_mddat(env,blk,mddat,istat) else if (blk%header == 'dynamics.meta') then - call parse_metadyn(blk,mddat) + call parse_metadyn(blk,mddat,istat) included = .true. end if end do @@ -1015,148 +952,105 @@ subroutine parse_dynamics_data(env,mddat,dict,included) end if return end subroutine parse_dynamics_data - subroutine parse_mddat(env,blk,mddat) + subroutine parse_mddat(env,blk,mddat,istat) implicit none type(systemdata),intent(inout) :: env type(datablock),intent(in) :: blk type(mddata),intent(inout) :: mddat + integer,intent(inout) :: istat logical,allocatable :: atlist(:) integer :: i,j,nat + logical :: rd if (blk%header .ne. 'dynamics') return nat = env%ref%nat allocate (atlist(nat),source=.false.) do i = 1,blk%nkv - call parse_md(mddat,blk%kv_list(i)) - - select case (blk%kv_list(i)%key) - case ('includermsd','atlist+') - call get_atlist(nat,atlist,blk%kv_list(i)%rawvalue,env%ref%at) - if (.not.allocated(env%includeRMSD)) allocate (env%includeRMSD(nat),source=1) - do j = 1,nat - if (atlist(j)) env%includeRMSD(j) = 1 - end do - - case ('excludermsd','atlist-') - call get_atlist(nat,atlist,blk%kv_list(i)%rawvalue,env%ref%at) - if (.not.allocated(env%includeRMSD)) allocate (env%includeRMSD(nat),source=1) - do j = 1,nat - if (atlist(j)) env%includeRMSD(j) = 0 - end do - - end select - + call parse_md(env,mddat,blk%kv_list(i),rd) + if (.not.rd) then + istat = istat+1 + write (stdout,fmturk) '['//blk%header//']-block',blk%kv_list(i)%key + end if end do deallocate (atlist) return end subroutine parse_mddat - subroutine parse_md_auto(mddat,kv) + subroutine parse_md_auto(env,mddat,kv,rd) implicit none + type(systemdata),intent(inout) :: env type(mddata) :: mddat type(keyvalue) :: kv - select case (kv%id) - case (1) !> float - call parse_md(mddat,kv%key,kv%value_f) - case (2) !> int - call parse_md(mddat,kv%key,kv%value_i) - case (3) !> bool - call parse_md(mddat,kv%key,kv%value_b) - case (4) !> string - call parse_md(mddat,kv%key,kv%value_c) - case default - select case (kv%key) - case ('active','active_levels') - mddat%active_potentials = kv%value_ia - end select - end select - end subroutine parse_md_auto - subroutine parse_md_float(mddat,key,val) - implicit none - type(mddata) :: mddat - character(len=*) :: key - real(wp) :: val - select case (key) + logical,intent(out) :: rd + logical,allocatable :: atlist(:) + integer :: nat,j + logical :: ex + rd = .true. + + select case (kv%key) + case ('active','active_levels') + mddat%active_potentials = kv%value_ia + + case ('includermsd','atlist+') + call get_atlist(nat,atlist,kv%rawvalue,env%ref%at) + if (.not.allocated(env%includeRMSD)) allocate (env%includeRMSD(nat),source=1) + do j = 1,nat + if (atlist(j)) env%includeRMSD(j) = 1 + end do + + case ('excludermsd','atlist-') + call get_atlist(nat,atlist,kv%rawvalue,env%ref%at) + if (.not.allocated(env%includeRMSD)) allocate (env%includeRMSD(nat),source=1) + do j = 1,nat + if (atlist(j)) env%includeRMSD(j) = 0 + end do + case ('length','length_ps') - mddat%length_ps = val + mddat%length_ps = kv%value_f case ('dump') - mddat%dumpstep = val + mddat%dumpstep = kv%value_f case ('hmass') - mddat%md_hmass = val + mddat%md_hmass = kv%value_f case ('tstep') - mddat%tstep = val + mddat%tstep = kv%value_f case ('t','temp','temperature') - mddat%tsoll = val + mddat%tsoll = kv%value_f mddat%thermostat = .true. - case default - return - end select - return - end subroutine parse_md_float - subroutine parse_md_int(mddat,key,val) - implicit none - type(mddata) :: mddat - character(len=*) :: key - integer :: val - real(wp) :: fval - select case (key) - case ('length','length_ps','dump','hmass','tstep') - fval = float(val) - call parse_md(mddat,key,fval) + case ('shake') - if (val <= 0) then - mddat%shake = .false. - else - mddat%shake = .true. - mddat%shk%shake_mode = min(val,2) - end if + select case (kv%id) + case (valuetypes%int) + if (kv%value_i <= 0) then + mddat%shake = .false. + else + mddat%shake = .true. + mddat%shk%shake_mode = min(kv%value_i,2) + end if + case (valuetypes%bool) + mddat%shake = kv%value_b + if (kv%value_b) mddat%shk%shake_mode = 1 + end select case ('printstep') - mddat%printstep = val - case('blocklength','blockl') - mddat%blockl = val - case ('t','temp','temperature') - mddat%tsoll = float(val) - mddat%thermostat = .true. - case default - return - end select - return - end subroutine parse_md_int - subroutine parse_md_c(mddat,key,val) - implicit none - type(mddata) :: mddat - character(len=*) :: key - character(len=*) :: val - logical :: ex - select case (key) - case('restart') - inquire(file=trim(val),exist=ex) - if(ex)then + mddat%printstep = kv%value_i + case ('blocklength','blockl') + mddat%blockl = kv%value_i + + case ('restart') + inquire (file=trim(kv%value_c),exist=ex) + if (ex) then mddat%restart = .true. - mddat%restartfile = trim(val) - endif - case default - return - end select - return - end subroutine parse_md_c - subroutine parse_md_bool(mddat,key,val) - implicit none - type(mddata) :: mddat - character(len=*) :: key - logical :: val - select case (key) - case ('shake') - mddat%shake = val - if (val) mddat%shk%shake_mode = 1 + mddat%restartfile = trim(kv%value_c) + end if + case default + rd = .false. return end select - return - end subroutine parse_md_bool + + end subroutine parse_md_auto !========================================================================================! - subroutine parse_metadyn(blk,mddat) + subroutine parse_metadyn(blk,mddat,istat) !************************************************** !* The following routines are used to !* read information into the "metadynamics" object @@ -1165,109 +1059,69 @@ subroutine parse_metadyn(blk,mddat) implicit none type(datablock),intent(in) :: blk type(mddata),intent(inout) :: mddat + integer,intent(inout) :: istat logical :: success type(mtdpot) :: mtd integer :: i,k + logical :: rd call mtd%deallocate() + success = .false. if (blk%header .ne. 'dynamics.meta') return do i = 1,blk%nkv - call parse_metadyn_auto(mtd,blk%kv_list(i),success) + call parse_metadyn_auto(mtd,blk%kv_list(i),success,rd) + if (.not.rd) then + istat = istat+1 + write (stdout,fmturk) '[['//blk%header//']]-block',blk%kv_list(i)%key + end if end do - call mddat%add(mtd) + if (success) call mddat%add(mtd) return end subroutine parse_metadyn - subroutine parse_metadyn_auto(mtd,kv,success) + subroutine parse_metadyn_auto(mtd,kv,success,rd) implicit none type(keyvalue) :: kv type(mtdpot) :: mtd - logical,intent(out) :: success - success = .false. - select case (kv%id) - case (1) !> float - call parse_mtd(mtd,kv%key,kv%value_f) - success = .true. - case (2) !> int - call parse_mtd(mtd,kv%key,kv%value_i) - success = .true. - case (3) !> bool - call parse_mtd(mtd,kv%key,kv%value_b) - success = .true. - case (4) !> string - call parse_mtd(mtd,kv%key,kv%value_c) - success = .true. - end select - end subroutine parse_metadyn_auto - subroutine parse_mtd_float(mtd,key,val) - implicit none - type(mtdpot) :: mtd - character(len=*) :: key - real(wp) :: val - select case (key) + logical,intent(inout) :: success + logical,intent(out) :: rd + rd = .true. + + select case (kv%key) +!>--- floats case ('alpha') - mtd%alpha = val + mtd%alpha = kv%value_f case ('kpush') - mtd%kpush = val + mtd%kpush = kv%value_f case ('dump','dump_fs') - mtd%cvdump_fs = val + mtd%cvdump_fs = kv%value_f case ('dump_ps') - mtd%cvdump_fs = val*1000.0_wp + mtd%cvdump_fs = kv%value_f*1000.0_wp case ('ramp') - mtd%ramp = val - case default - return - end select - return - end subroutine parse_mtd_float - subroutine parse_mtd_int(mtd,key,val) - implicit none - type(mtdpot) :: mtd - character(len=*) :: key - integer :: val - real(wp) :: fval - select case (key) - case ('type') - mtd%mtdtype = val - case ('dump','dump_fs','dump_ps') - fval = float(val) - call parse_mtd(mtd,key,fval) - case default - return - end select - return - end subroutine parse_mtd_int - subroutine parse_mtd_c(mtd,key,val) - implicit none - type(mtdpot) :: mtd - character(len=*) :: key - character(len=*) :: val - select case (key) + mtd%ramp = kv%value_f + +!>--- strings case ('type') - select case (val) - case ('rmsd') - mtd%mtdtype = cv_rmsd - case default - mtd%mtdtype = 0 + success = .true. + select case (kv%id) + case (valuetypes%int) + mtd%mtdtype = kv%value_i + case (valuetypes%string) + select case (kv%value_c) + case ('rmsd') + mtd%mtdtype = cv_rmsd + case default + mtd%mtdtype = 0 + end select end select case ('biasfile') mtd%mtdtype = cv_rmsd_static - mtd%biasfile = val + mtd%biasfile = kv%value_c case default + rd = .false. return end select - return - end subroutine parse_mtd_c - subroutine parse_mtd_bool(mtd,key,val) - implicit none - type(mtdpot) :: mtd - character(len=*) :: key - logical :: val - select case (key) - case default - return - end select - return - end subroutine parse_mtd_bool + + end subroutine parse_metadyn_auto !========================================================================================! !========================================================================================! diff --git a/src/parsing/parse_datastruct.f90 b/src/parsing/parse_datastruct.f90 index 461b87e7..7dd7e95d 100644 --- a/src/parsing/parse_datastruct.f90 +++ b/src/parsing/parse_datastruct.f90 @@ -44,6 +44,7 @@ module parse_datastruct procedure :: addblk => root_addblk procedure :: lowercase_keys => root_lowercase_keys procedure :: print => root_print + procedure :: print2 => root_print2 end type root_object !========================================================================================! @@ -124,5 +125,24 @@ subroutine root_print(self) end do end subroutine root_print + + subroutine root_print2(self) + class(root_object) :: self + integer :: i + if (allocated(self%filename)) then + write(stdout,'(a)') repeat('*',80) + write (stdout,'("*",1x,a,a,a)') 'INPUT FILE ',self%filename,' content (without comments):' + end if + write(stdout,'(a)') repeat('*',80) + do i = 1,self%nkv + write (stdout,'("*",1x,a)') trim(self%kv_list(i)%print2()) + end do + do i = 1,self%nblk + call self%blk_list(i)%print2() + end do + write(stdout,'(a)') repeat('*',80) + end subroutine root_print2 + + !========================================================================================! end module parse_datastruct diff --git a/src/parsing/parse_keyvalue.f90 b/src/parsing/parse_keyvalue.f90 index e33a888e..677b7686 100644 --- a/src/parsing/parse_keyvalue.f90 +++ b/src/parsing/parse_keyvalue.f90 @@ -40,6 +40,7 @@ module parse_keyvalue character(len=:),allocatable :: value_ca(:) !> id=8, an array of strings/multiline string contains procedure :: print => print_kv + procedure :: print2 => print_kv2 procedure :: deallocate => deallocate_kv procedure :: set_valuestring => kv_set_valuestring procedure :: add_raw_array_string => kv_add_raw_array_string @@ -374,17 +375,20 @@ subroutine get_valueautotype(kv) return end if - !> real if (is_float(atmp)) then + !>--- real read (atmp,*,iostat=io) num kv%id = valuetypes%float kv%value_f = num + kv%value_i = nint(num) !> backed up - !> integer else if (is_int(atmp)) then + !>--- integer read (atmp,*,iostat=io) num kv%id = valuetypes%int kv%value_i = nint(num) + kv%value_f = real(num,wp) !> backed up + end if return @@ -443,13 +447,14 @@ subroutine kv_set_valuestring(self) class(keyvalue) :: self character(len=20) :: atmp character(len=:),allocatable :: btmp + integer :: i if (allocated(self%rawvalue)) deallocate (self%rawvalue) select case (self%id) case (valuetypes%float) !> float - write (atmp,'(e20.10)') self%value_f + write (atmp,'(f20.10)') self%value_f atmp = adjustl(atmp) case (valuetypes%int) !> integer - write (atmp,'(i20)') self%value_i + write (atmp,'(i0)') self%value_i atmp = adjustl(atmp) case (valuetypes%bool) !> boolean if (self%value_b) then @@ -461,26 +466,53 @@ subroutine kv_set_valuestring(self) btmp = self%value_c case (valuetypes%float_array) !> float array - write (atmp,'(e20.5)') self%value_fa(1) - btmp = '['//trim(adjustl(atmp))//', ... ]' + btmp = '[' + do i = 1,self%na-1 + write (atmp,*) self%value_fa(i) + btmp = trim(btmp)//trim(atmp)//',' + end do + write (atmp,*) self%value_fa(self%na) + btmp = trim(btmp)//trim(atmp)//']' case (valuetypes%int_array) !> int array - write (atmp,'(i20)') self%value_ia(1) - btmp = '['//trim(adjustl(atmp))//', ... ]' + btmp = '[' + do i = 1,self%na-1 + write (atmp,'(i0)') self%value_ia(i) + btmp = trim(btmp)//trim(atmp)//',' + end do + write (atmp,'(i0)') self%value_ia(self%na) + btmp = trim(btmp)//trim(atmp)//']' case (valuetypes%bool_array) !> bool array - if (self%value_ba(1)) then + btmp = '[' + do i = 1,self%na-1 + if (self%value_ba(i)) then + atmp = 'true' + else + atmp = 'false' + end if + btmp = trim(btmp)//trim(atmp)//',' + end do + if (self%value_ba(self%na)) then atmp = 'true' else atmp = 'false' end if - btmp = '['//trim(adjustl(atmp))//', ... ]' + btmp = trim(btmp)//trim(atmp)//']' case (valuetypes%string_array) !> multiline comment - btmp = '["'//trim(self%value_ca(1))//'", ..., etc.' + btmp = '[' + do i = 1,self%na-1 + btmp = trim(btmp)//'"'//trim(self%value_rawa(i))//'",' + end do + btmp = trim(btmp)//'"'//trim(self%value_rawa(self%Na))//'"]' case (valuetypes%raw_array) !> unspecified array - btmp = '['//trim(self%value_rawa(1))//', ..., etc.' + btmp = '[' + do i = 1,self%na-1 + btmp = trim(btmp)//trim(self%value_rawa(i))//',' + end do + btmp = trim(btmp)//trim(self%value_rawa(self%Na))//']' case default atmp = '' @@ -611,5 +643,18 @@ subroutine truncate15(ins,outs) outs = ins(1:l) end if end subroutine truncate15 + + function print_kv2(self) result(outstr) + implicit none + character(len=:),allocatable :: outstr + class(keyvalue) :: self + select case (self%id) + case (valuetypes%string) + outstr = self%key//' = "'//self%rawvalue//'"' + case default + outstr = self%key//' = '//self%rawvalue + end select + end function print_kv2 + !========================================================================================! end module parse_keyvalue diff --git a/src/parsing/parse_maindata.f90 b/src/parsing/parse_maindata.f90 index a3b90e56..723ddc1f 100644 --- a/src/parsing/parse_maindata.f90 +++ b/src/parsing/parse_maindata.f90 @@ -37,37 +37,56 @@ module parse_maindata implicit none public + character(len=*),parameter,private :: fmturk = '("unrecognized KEYWORD in ",a," : ",a)' + character(len=*),parameter,private :: fmtura = '("unrecognized ARGUMENT : ",a)' + !========================================================================================! !========================================================================================! contains !> MODULE PROCEDURES START HERE !========================================================================================! !========================================================================================! - subroutine parse_main_auto(env,kv) + subroutine parse_main_auto(env,kv,istat) implicit none type(systemdata) :: env type(keyvalue) :: kv + integer,intent(inout) :: istat + integer :: istat_ref + logical :: rd + istat_ref = istat + rd = .false. select case (kv%id) - case (valuetypes%float) !> float - call parse_main_float(env,kv%key,kv%value_f) - case (valuetypes%int) !> int - call parse_main_int(env,kv%key,kv%value_i) - case (valuetypes%bool) !> bool - call parse_main_bool(env,kv%key,kv%value_b) - case (valuetypes%string) !> string - call parse_main_c(env,kv%key,kv%value_c) - end select -!> other, with multiple or raw type - select case (kv%key) - case ('optlev','ancopt_level') - env%optlev = optlevnum(kv%rawvalue) + case (valuetypes%float) !>--- float + call parse_main_float(env,kv%key,kv%value_f,rd) + case (valuetypes%int) !>--- int + call parse_main_int(env,kv%key,kv%value_i,rd) + case (valuetypes%bool) !>--- bool + call parse_main_bool(env,kv%key,kv%value_b,rd) + case (valuetypes%string) !>--- string + call parse_main_c(env,kv%key,kv%value_c,rd) end select +!>--- other, with multiple or raw type + if (.not.rd) then + select case (kv%key) + case ('optlev','ancopt_level') + env%optlev = optlevnum(kv%rawvalue) + case default + istat = istat+1 + end select + end if +!>--- if none of the options was recognizeda and istat increased as a consequence, print that + if (istat > istat_ref) then + write (stdout,fmturk) 'main section',kv%key + end if + end subroutine parse_main_auto - subroutine parse_main_float(env,key,val) + subroutine parse_main_float(env,key,val,rd) implicit none type(systemdata) :: env character(len=*) :: key real(wp) :: val + logical,intent(out) :: rd + rd = .true. select case (key) case ('wscal') env%potscal = val @@ -75,29 +94,36 @@ subroutine parse_main_float(env,key,val) case ('wpad') env%potpad = val env%wallsetup = .true. - + case default + rd = .false. end select return end subroutine parse_main_float - subroutine parse_main_int(env,key,val) + subroutine parse_main_int(env,key,val,rd) implicit none type(systemdata) :: env character(len=*) :: key + logical,intent(out) :: rd integer :: val + rd = .true. select case (key) case ('threads','parallel') env%Threads = val env%autothreads = .true. env%threadssetmanual = .true. + case default + rd = .false. end select return end subroutine parse_main_int - subroutine parse_main_c(env,key,val) + subroutine parse_main_c(env,key,val,rd) implicit none type(systemdata) :: env character(len=*) :: key character(len=*) :: val + logical,intent(out) :: rd type(coord) :: mol + rd = .true. select case (key) case ('bin','binary') env%ProgName = val @@ -119,9 +145,9 @@ subroutine parse_main_c(env,key,val) env%preopt = .false. env%crestver = crest_optimize env%optlev = 0.0_wp - if(val.eq.'ohess')then - env%crest_ohess=.true. - endif + if (val .eq. 'ohess') then + env%crest_ohess = .true. + end if case ('ancopt_ensemble','optimize_ensemble','mdopt') env%preopt = .false. env%crestver = crest_mdopt2 @@ -178,8 +204,16 @@ subroutine parse_main_c(env,key,val) env%preopt = .false. env%crestver = crest_rigcon env%runver = crest_rigcon + + case ('protonate') + env%properties = p_protonate + env%crestver = crest_protonate + case default - env%crestver = crest_imtd + !>--- keyword was recognized, but invalid argument supplied + write (stdout,fmtura) val + error stop + end select case ('ensemble_input','ensemble','input_ensemble') env%ensemblename = val @@ -197,14 +231,18 @@ subroutine parse_main_c(env,key,val) case ('watlist','wat') env%potatlist = val env%wallsetup = .true. + case default + rd = .false. end select return end subroutine parse_main_c - subroutine parse_main_bool(env,key,val) + subroutine parse_main_bool(env,key,val,rd) implicit none type(systemdata) :: env character(len=*) :: key + logical,intent(out) :: rd logical :: val + rd = .true. select case (key) case ('preopt') env%preopt = val @@ -219,18 +257,20 @@ subroutine parse_main_bool(env,key,val) call read_restart(env) end if case ('multilevelopt') - env%multilevelopt = val + env%multilevelopt = val case ('refine_presort') env%refine_presort = val case ('omp_nested') env%omp_allow_nested = val + case default + rd = .false. end select return end subroutine parse_main_bool !========================================================================================! - subroutine parse_main_blk(env,blk) + subroutine parse_main_blk(env,blk,istat) !************************************** !* Some shorter blocks are not defined !* in separate source files. They can @@ -239,18 +279,21 @@ subroutine parse_main_blk(env,blk) implicit none type(systemdata) :: env type(datablock) :: blk + integer,intent(inout) :: istat select case (blk%header) case ('cregen') - call parse_cregen(env,blk) + call parse_cregen(env,blk,istat) case ('confsolv') - call parse_confsolv(env,blk) + call parse_confsolv(env,blk,istat) case ('thermo') - call parse_thermo(env,blk) + call parse_thermo(env,blk,istat) + case ('protonation') + call parse_protonation(env,blk,istat) end select end subroutine parse_main_blk !========================================================================================! - subroutine parse_cregen(env,blk) + subroutine parse_cregen(env,blk,istat) !**************************************** !* parse settings for the CREGEN routine !**************************************** @@ -258,6 +301,7 @@ subroutine parse_cregen(env,blk) type(systemdata) :: env type(datablock) :: blk type(keyvalue) :: kv + integer,intent(inout) :: istat integer :: i !>--- parse the arguments do i = 1,blk%nkv @@ -273,21 +317,25 @@ subroutine parse_cregen(env,blk) env%bthr2 = kv%value_f case ('eqv','nmr') env%doNMR = kv%value_b + case default + !>--- unrecognized keyword + istat = istat+1 + write (stdout,fmturk) '[cregen]-block',kv%key end select end do end subroutine parse_cregen !========================================================================================! - subroutine parse_confsolv(env,blk) + subroutine parse_confsolv(env,blk,istat) use ConfSolv_module implicit none type(systemdata) :: env type(datablock) :: blk type(keyvalue) :: kv + integer,intent(inout) :: istat integer :: i !>--- add ConfSolv as refinement level to give a ΔΔGsoln call env%addrefine(refine%ConfSolv) - !env%refine_presort = .true. env%ewin = 100.0_wp !>--- parse the arguments @@ -307,12 +355,12 @@ subroutine parse_confsolv(env,blk) if (kv%na == 2) then cs_solvent = trim(kv%value_rawa(1)) cs_smiles = trim(kv%value_rawa(2)) - else if(index(kv%value_c,'.csv').ne.0)then + else if (index(kv%value_c,'.csv') .ne. 0) then cs_solvfile = kv%value_c else cs_solvent = kv%value_c end if - case('solvent_csv','solvfile') + case ('solvent_csv','solvfile') cs_solvfile = kv%value_c case ('solvent_name') cs_solvent = kv%value_c @@ -320,13 +368,17 @@ subroutine parse_confsolv(env,blk) cs_smiles = kv%value_c case ('model_path','param','checkpoint') cs_param = kv%value_c + case default + !>--- unrecognized keyword + istat = istat+1 + write (stdout,fmturk) '[confsolv]-block',kv%key end select end do end subroutine parse_confsolv !========================================================================================! - subroutine parse_thermo(env,blk) + subroutine parse_thermo(env,blk,istat) !**************************************** !* parse settings for the Thermo routine !**************************************** @@ -334,6 +386,7 @@ subroutine parse_thermo(env,blk) type(systemdata) :: env type(datablock) :: blk type(keyvalue) :: kv + integer,intent(inout) :: istat integer :: i !>--- parse the arguments do i = 1,blk%nkv @@ -348,24 +401,71 @@ subroutine parse_thermo(env,blk) case ('trange') if (kv%na >= 2) then env%thermo%trange(1) = minval(kv%value_fa(1:2),1) - env%thermo%trange(2) = maxval(kv%value_fa(1:2),1) - endif - if(kv%na >= 3)then + env%thermo%trange(2) = maxval(kv%value_fa(1:2),1) + end if + if (kv%na >= 3) then env%thermo%trange(3) = kv%value_fa(3) - endif + end if case ('tstep') env%thermo%trange(3) = kv%value_f case ('input','coords') env%thermo%coords = kv%value_c - if(allocated(env%thermo%vibfile)) env%properties = p_thermo + if (allocated(env%thermo%vibfile)) env%properties = p_thermo case ('freq_input','vibs','hessian') - env%thermo%vibfile = kv%value_c - if(allocated(env%thermo%coords)) env%properties = p_thermo + env%thermo%vibfile = kv%value_c + if (allocated(env%thermo%coords)) env%properties = p_thermo + case default + !>--- unrecognized keyword + istat = istat+1 + write (stdout,fmturk) '[thermo]-block',kv%key end select end do end subroutine parse_thermo +!========================================================================================! + subroutine parse_protonation(env,blk,istat) +!****************************************** +!* parse settings for protonation settings +!****************************************** + implicit none + type(systemdata) :: env + type(datablock) :: blk + type(keyvalue) :: kv + integer,intent(inout) :: istat + integer :: i + external :: swparse +!>--- parse the arguments + do i = 1,blk%nkv + kv = blk%kv_list(i) + select case (kv%key) + case ('ewin') + env%protb%ewin = kv%value_f + case ('swel','ion') + call swparse(kv%value_c,env%protb) + case ('ffopt') + env%protb%ffopt = kv%value_b + case ('freezeopt') + env%protb%hnewopt = kv%value_b + case ('finalopt') + env%protb%finalopt = kv%value_b + + case ('activelmo') + env%protb%active_lmo(1:) = kv%value_ba(1:) + case ('pi') + env%protb%active_lmo(1) = kv%value_b + case ('lp') + env%protb%active_lmo(2) = kv%value_b + case ('delpi','delocpi') + env%protb%active_lmo(3) = kv%value_b + case default + !>--- unrecognized keyword + istat = istat+1 + write (stdout,fmturk) '[protonation]-block',kv%key + end select + end do + end subroutine parse_protonation + !========================================================================================! !========================================================================================! end module parse_maindata diff --git a/src/qmhelpers/local.f90 b/src/qmhelpers/local.f90 index eaca0953..26915cbc 100644 --- a/src/qmhelpers/local.f90 +++ b/src/qmhelpers/local.f90 @@ -103,6 +103,9 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z, & logical l1,l2,l3,flip integer LWORK,IWORK,LIWORK,INFO integer :: iscreen,icoord,ilmoi,icent ! file handles + integer :: ntmpsave + real(wp),allocatable :: tmpsave(:,:) + !> BLAS external dgemm @@ -275,6 +278,8 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z, & call mocent(nat,nao2,n,cmo,s,qmo,xcen,aoat2) allocate (rklmo(5,2*n)) + allocate (tmpsave(4,2*n), source=0.0_wp) + ntmpsave=0 !if (pr_local) write (*,*) 'lmo centers(Z=2) and atoms on file ' if (pr_local) write (*,*) ' LMO type Fii/eV ncent charge center contributions...' ! if (pr_local) open (newunit=iscreen,file='xtbscreen.xyz') @@ -330,7 +335,7 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z, & !>--- 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 +! 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 @@ -346,10 +351,17 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z, & ! 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 if +!>--- tmpsave + if(jdum.gt.1)then + ntmpsave = ntmpsave + 1 + tmpsave(1:3, ntmpsave) = ecent(i,1:3) + tmpsave(4,ntmpsave) = float(jdum) + endif + end do !ccccccccccccccccccccccccccccccccccccccccccccccccccccccc @@ -471,6 +483,11 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z, & ! end do ! write (iscreen,'(a2,3F24.10)') i2e(1),dtot(1:3)*autoaa end if + !>--- tmpsave + ntmpsave = ntmpsave + 1 + tmpsave(1:3, ntmpsave) = dtot(1:3) + tmpsave(4,ntmpsave) = rklmo(5,pmo) + imo = lneigh(1,pmo) vec1(1:3) = xyz(1:3,nm) @@ -482,7 +499,6 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z, & 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,*) @@ -491,6 +507,10 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z, & ! end do ! write (iscreen,'(a2,3F24.10)') i2e(1),dtot(1:3)*autoaa end if + !>--- tmpsave + ntmpsave = ntmpsave + 1 + tmpsave(1:3, ntmpsave) = dtot(1:3) + tmpsave(4,ntmpsave) = rklmo(5,imo) m = m+1 end do @@ -506,22 +526,21 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z, & !ccccccccccccccccccccccccccccccccccccccccccccccccccccccc !> collect possible protonation sites for output - nprot = 0 - k = new+n - do i=1,k - if(rklmo(5,i) > 1) nprot=nprot+1 + nprot = ntmpsave + do i =1,ntmpsave + tmpsave(4,i) = tmpsave(4,i) - 1.0_wp !> adjust numbering (we discard sigma anyways) + if(tmpsave(4,i) < 1.0_wp) nprot=nprot-1 enddo if(nprot > 0)then allocate(protxyz(4,nprot), source=0.0_wp) - m=0 - do i=1,k - 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 + ii=0 + do i=1,ntmpsave + if(tmpsave(4,i) < 1.0_wp) cycle + ii=ii+1 + protxyz(:,ii) = tmpsave(:,i) + enddo endif + deallocate(tmpsave) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccc @@ -556,14 +575,6 @@ subroutine local(nat,at,nbf,nao,ihomo,xyz,z, & ! 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 diff --git a/src/sortens.f90 b/src/sortens.f90 index bcc2cc3c..c377571f 100644 --- a/src/sortens.f90 +++ b/src/sortens.f90 @@ -54,7 +54,7 @@ subroutine sort_ens(sort,infile,verbose) !sort=env%protb - associate (popthr => sort%popthr, & + associate (& & allowFrag => sort%allowFrag,ewin => sort%ewin, & & protdeprot => sort%protdeprot,deprotprot => sort%deprotprot, & & nfrag => sort%nfrag,threshsort => sort%threshsort) diff --git a/src/strucreader.f90 b/src/strucreader.f90 index 56078d52..b2272e71 100644 --- a/src/strucreader.f90 +++ b/src/strucreader.f90 @@ -263,7 +263,7 @@ subroutine rdensembleparam(fname,nat,nall,conform) character(len=*),intent(in) :: fname integer,intent(out) :: nat integer,intent(out) :: nall - logical,optional :: conform + logical,intent(out),optional :: conform logical :: conformdum integer :: dum,iosum integer :: natref @@ -535,7 +535,7 @@ subroutine rdensemble_mixed2(fname,natmax,nall,nats,ats,xyz) end subroutine rdensemble_mixed2 !========================================================================================! - subroutine rdensemble_coord_type(fname,nall,ensemble) + subroutine rdensemble_coord_type(fname,nall,structures) !********************************************************* !* subroutine rdensemble_coord_type !* A variant of the rdensemble routine that automatically @@ -544,34 +544,59 @@ subroutine rdensemble_coord_type(fname,nall,ensemble) implicit none character(len=*),intent(in) :: fname !> name of the ensemble file integer,intent(out) :: nall !> number of structures in ensemble - type(coord),intent(out),allocatable :: ensemble(:) + type(coord),intent(out),allocatable :: structures(:) real(wp),allocatable :: xyz(:,:,:) integer :: nat + integer,allocatable :: nats(:) integer,allocatable :: at(:) + integer,allocatable :: ats(:,:) real(wp),allocatable :: eread(:) - integer :: i,j,k,ich,io - logical :: ex - - call rdensembleparam(fname,nat,nall) - allocate (xyz(3,nat,nall),at(nat),eread(nall)) - call rdensemble(fname,nat,nall,at,xyz,eread) - !>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<--- Important: coord types must be in Bohrs - xyz = xyz/bohr - !>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<--- all the same molecule +! allocate (xyz(3,nat,nall),at(nat),eread(nall)) +! call rdensemble(fname,nat,nall,at,xyz,eread) +! !>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<--- Important: coord types must be in Bohrs +! xyz = xyz/bohr +! !>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<--- multiple sizes + allocate(structures(nall)) + allocate(xyz(3,nat,nall),ats(nat,nall),nats(nall),eread(nall)) + call rdensemble_mixed2(fname,nat,nall,nats,ats,xyz) + !>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<--- Important: coord types must be in Bohrs + xyz = xyz/bohr + !>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<