diff --git a/src/calculator/tblite_api.F90 b/src/calculator/tblite_api.F90 index 0790aac5..3e7a4edf 100644 --- a/src/calculator/tblite_api.F90 +++ b/src/calculator/tblite_api.F90 @@ -36,7 +36,7 @@ module tblite_api use tblite_param, only : param_record use tblite_results,only:tblite_resultstype => results_type use tblite_wavefunction_mulliken,only:get_molecular_dipole_moment - use tblite_ceh_singlepoint,only:ceh_guess + use tblite_ceh_singlepoint,only:ceh_singlepoint use tblite_ceh_ceh,only:new_ceh_calculator #endif use wiberg_mayer @@ -137,19 +137,19 @@ subroutine tblite_setup(mol,chrg,uhf,lvl,etemp,tblite) select case (tblite%lvl) case (xtblvl%gfn1) if (pr) call tblite%ctx%message("tblite> setting up GFN1-xTB calculation") - call new_gfn1_calculator(tblite%calc,mctcmol) + call new_gfn1_calculator(tblite%calc,mctcmol,error) case (xtblvl%gfn2) if (pr) call tblite%ctx%message("tblite> setting up GFN2-xTB calculation") - call new_gfn2_calculator(tblite%calc,mctcmol) + call new_gfn2_calculator(tblite%calc,mctcmol,error) case (xtblvl%ipea1) if (pr) call tblite%ctx%message("tblite> setting up IPEA1-xTB calculation") - call new_ipea1_calculator(tblite%calc,mctcmol) + call new_ipea1_calculator(tblite%calc,mctcmol,error) case (xtblvl%ceh) if (pr) call tblite%ctx%message("tblite> setting up CEH calculation") - call new_ceh_calculator(tblite%calc,mctcmol) + call new_ceh_calculator(tblite%calc,mctcmol,error) case (xtblvl%eeq) if (pr) call tblite%ctx%message("tblite> setting up D4 EEQ charges calculation") - call new_ceh_calculator(tblite%calc,mctcmol) !> doesn't matter but needs initialization + call new_ceh_calculator(tblite%calc,mctcmol,error) !> doesn't matter but needs initialization case (xtblvl%param) if (pr) call tblite%ctx%message("tblite> setting up xtb calculator from parameter file") if(allocated(tblite%paramfile))then @@ -161,7 +161,7 @@ subroutine tblite_setup(mol,chrg,uhf,lvl,etemp,tblite) endif else if (pr) call tblite%ctx%message("tblite> parameter file does not exist, defaulting to GFN2-xTB") - call new_gfn2_calculator(tblite%calc,mctcmol) + call new_gfn2_calculator(tblite%calc,mctcmol,error) endif case default call tblite%ctx%message("Error: Unknown method in tblite!") @@ -318,7 +318,7 @@ subroutine tblite_singlepoint(mol,chrg,uhf,tblite,energy,gradient,iostatus) & energy,gradient, & & sigma,verbosity,results=tblite%res) case (xtblvl%ceh) - call ceh_guess(tblite%ctx,tblite%calc,mctcmol,error,tblite%wfn, & + call ceh_singlepoint(tblite%ctx,tblite%calc,mctcmol,tblite%wfn, & & tblite%accuracy,verbosity) case(xtblvl%eeq) call eeq_guess(mctcmol, tblite%calc, tblite%wfn) diff --git a/src/cleanup.f90 b/src/cleanup.f90 index df4b1468..d9071e0e 100644 --- a/src/cleanup.f90 +++ b/src/cleanup.f90 @@ -52,6 +52,34 @@ end subroutine rmrfw !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC !C Specific cleanup routines for different parts of the CREST code !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC + +subroutine custom_cleanup(env) + use crest_data + use iomod + implicit none + type(systemdata) :: env + integer :: i + if(.not.env%keepModef)then + call rmrfw('METADYN') + call rmrfw('NORMMD') + call rmrf('MRMSD') + call rmrf('TRIALMD') + call rmrf('TRIALOPT') + call rmrf('MDFILES') + call rmrf('OPTIM') + call rmrf('PROP') + call rmrfw('.cre_') + call rmrf('cregen_*.tmp') + call rmrf('MDFILES') + if(.not.any(env%calc%calcs(:)%pr))then + call rmrfw('calculation.level.') + endif + endif + call rmrf('.CHRG .UHF') + call rmrf('.history.xyz') +end subroutine custom_cleanup + + !------------------------------------------------------------------------- ! General cleanup function. Wipes most files that may be written by crest !------------------------------------------------------------------------- diff --git a/src/confparse.f90 b/src/confparse.f90 index 19bd094a..74a60a3b 100644 --- a/src/confparse.f90 +++ b/src/confparse.f90 @@ -250,20 +250,19 @@ subroutine parseflags(env,arg,nra) inquire (file='.CHRG',exist=ex) if (any(index(arg,'-chrg') .ne. 0)) ex = .false. if (ex) then - call rdshort('.CHRG',env%chrg) - if (env%chrg .ne. 0) then - write (*,'(2x,a,i0)') 'molecular charge read from .CHRG: ',env%chrg - end if + write(stdout,*) '**ERROR** CREST will not read .CHRG files following version 3.0.2' + write(stdout,*) 'Please use --chrg or the input file specifications.' + error stop end if inquire (file='.UHF',exist=ex) if (any(index(arg,'-uhf') .ne. 0)) ex = .false. if (ex) then - call rdshort('.UHF',env%uhf) - if (env%uhf .ne. 0) then - write (*,'(2x,a,i0)') 'nα-nβ electrons read from .UHF: ',env%uhf - end if + write(stdout,*) '**ERROR** CREST will not read .UHF files following version 3.0.2' + write(stdout,*) 'Please use --uhf or the input file specifications.' + error stop end if + !>--- options for constrained conformer sampling env%fixfile = 'none selected' diff --git a/src/cregen.f90 b/src/cregen.f90 index 6004d5c5..905cfb1c 100644 --- a/src/cregen.f90 +++ b/src/cregen.f90 @@ -78,6 +78,8 @@ subroutine newcregen(env,quickset,infile) character(len=128),allocatable :: comments(:) character(len=128),allocatable :: comref(:) real(wp),allocatable :: er(:) !> energies + type(coord),allocatable :: structures(:) !> a list of structures using the coord type + type(coord),allocatable :: references(:) !> the reference structure list !>--- dummy ensemble arguments integer :: nallref integer :: nallnew @@ -159,6 +161,9 @@ subroutine newcregen(env,quickset,infile) !>--- allocate space and read in the ensemble allocate (at(nat),comments(nallref),xyz(3,nat,nallref)) call rdensemble(fname,nat,nallref,at,xyz,comments) + !call rdensemble(fname,nallref,structures) + !allocate(references, source=structures) + !>--- track ensemble for restart call trackensemble(fname,nat,nallref,at,xyz,comments) @@ -215,7 +220,7 @@ subroutine newcregen(env,quickset,infile) !>--- do the rotational constants and RMSD check if (sortRMSD) then allocate (group(0:nall)) - call cregen_CRE(prch,env,nat,nall,at,xyz,comments,nallnew,group) + call cregen_CRE(prch,env,nat,nall,at,xyz,comments,nallnew,group,.false.) !>--- if structures were discarded, resize xyz if (nallnew .lt. nall) then nall = nallnew @@ -239,7 +244,7 @@ subroutine newcregen(env,quickset,infile) end if if (sortRMSD2) then allocate (group(0:nall)) - call cregen_CRE_2(prch,env,nat,nall,at,xyz,comments,nallnew,group) + call cregen_CRE(prch,env,nat,nall,at,xyz,comments,nallnew,group,.true.) end if !=====================================================================! @@ -413,7 +418,7 @@ subroutine cregen_prout(env,simpleset,pr1,pr2,pr3,pr4) pr3 = .false. !> plain energy list pr4 = .false. !> group list printout - if (any(simpleset == (/6,7/))) then + if (any(simpleset == (/6,7/)).or.env%esort) then pr1 = .false. pr2 = .false. if (env%crestver .ne. crest_solv) pr3 = .true. @@ -489,7 +494,7 @@ subroutine cregen_director(env,simpleset,checkbroken,sortE,sortRMSD,sortRMSD2, & anal = .false. end if - if (any(simpleset == (/6,7/))) then !energy sorting only + if (any(simpleset == (/6,7/)).or.env%esort) then !energy sorting only checkbroken = .false. sortE = .true. sortRMSD = .false. @@ -528,7 +533,7 @@ subroutine cregen_director(env,simpleset,checkbroken,sortE,sortRMSD,sortRMSD2, & anal = .false. end if - if(simpleset == 13)then !msreact mode + if (simpleset == 13) then !msreact mode checkbroken = .false. sorte = .true. sortRMSD = .true. @@ -538,10 +543,10 @@ subroutine cregen_director(env,simpleset,checkbroken,sortE,sortRMSD,sortRMSD2, & conffile = .true. topocheck = .false. checkez = .false. - bonusfiles= .false. + bonusfiles = .false. anal = .false. - endif + end if return end subroutine cregen_director @@ -1113,7 +1118,7 @@ end subroutine cregen_esort !=========================================================================================! -subroutine cregen_CRE(ch,env,nat,nall,at,xyz,comments,nallout,group) +subroutine cregen_CRE(ch,env,nat,nall,at,xyz,comments,nallout,group,nosort) !************************************************************* !* subroutine cregen_CRE !* sort the ensemble based on rotational constants,RMSD and @@ -1124,6 +1129,7 @@ subroutine cregen_CRE(ch,env,nat,nall,at,xyz,comments,nallout,group) !* at - atom types !* xyz - Cartesian coordinates !* comments - commentary lines containing the energy +!* nosort - don't actually sort !* On Output: nallout - new total number of structures !* group - to which group every structure belongs !************************************************************** @@ -1141,6 +1147,7 @@ subroutine cregen_CRE(ch,env,nat,nall,at,xyz,comments,nallout,group) integer,intent(in) :: at(nat) real(wp),intent(inout) :: xyz(3,nat,nall) integer,intent(inout) :: group(0:nall) + logical,intent(in) :: nosort character(len=*) :: comments(nall) integer :: nallout logical :: enantio = .true. !check for enantiomers? @@ -1377,29 +1384,35 @@ subroutine cregen_CRE(ch,env,nat,nall,at,xyz,comments,nallout,group) nallout = nall-k deallocate (mask) write (ch,*) 'number of doubles removed by rot/RMSD :',k - write (ch,*) 'total number unique points considered further :',nallout + if (.not.nosort) then + write (ch,*) 'total number unique points considered further :',nallout + else + nallout = nall + end if !>-- sort structures, energies, rot const. and comments - j = 0 - l = nall+1 - do i = 1,nall - if (double(i) .eq. 0) then - j = j+1 - orderref(i) = j - else - l = l-1 - orderref(i) = l - end if - end do - order = orderref - allocate (cdum(3,nat)) - call xyzqsort(nat,nall,xyz,cdum,order,1,nall) - deallocate (cdum) - order = orderref - call maskqsort(er,1,nall,order) - order = orderref - call stringqsort(nall,comments,1,nall,order) - order = orderref - call matqsort(3,nall,rot,rotdum,1,nall,order) + if (.not.nosort) then + j = 0 + l = nall+1 + do i = 1,nall + if (double(i) .eq. 0) then + j = j+1 + orderref(i) = j + else + l = l-1 + orderref(i) = l + end if + end do + order = orderref + allocate (cdum(3,nat)) + call xyzqsort(nat,nall,xyz,cdum,order,1,nall) + deallocate (cdum) + order = orderref + call maskqsort(er,1,nall,order) + order = orderref + call stringqsort(nall,comments,1,nall,order) + order = orderref + call matqsort(3,nall,rot,rotdum,1,nall,order) + end if !>-- finally, determine conformer groups and their rotamers allocate (c1(3,nat)) @@ -1450,349 +1463,11 @@ subroutine cregen_CRE(ch,env,nat,nall,at,xyz,comments,nallout,group) end if end do group(0) = k !>-- total number of groups - - deallocate (enuc,c1,double) - deallocate (order,orderref) - if (allocated(ecoul)) deallocate (ecoul) - deallocate (er) - deallocate (rot) - return -contains - function counth(nat,at) result(nh) - implicit none - integer :: nat - integer :: at(nat) - integer :: nh,i - nh = 0 - do i = 1,nat - if (at(i) == 1) nh = nh+1 - end do - return - end function counth - subroutine heavymask(nat,at,mask) - implicit none - integer :: nat - integer :: at(nat) - integer :: mask(nat) - integer :: i - mask = 0 - do i = 1,nat - if (at(i) .ne. 1) mask(i) = 1 - end do - return - end subroutine heavymask -end subroutine cregen_CRE - -!=========================================================================================! - -subroutine calc_ecoul(nat,at,xyz,cn,atommask,ecoulomb) -!******************************************************* -!* subroutine calc_ecoul -!* calculate a (artificial) coulomb energy for a given -!* structure. Some atoms can be ignored via a mask -!******************************************************* - use crest_parameters - implicit none - integer :: nat - integer :: at(nat) - real(wp) :: xyz(3,nat) !should be in Bohrs - real(wp) :: cn(nat) - integer :: atommask(nat) !contains 1 or 0 - real(wp) :: ecoulomb - integer :: i,j - real(wp) :: cnexp,rexp,natexp - real(wp) :: r,dum,zi,zj - ecoulomb = 0.0_wp - !cn=1.0_wp !comment in to ignore cn - cnexp = 1.0_wp - rexp = 1.0_wp - natexp = 3.5_wp - do i = 1,nat - if (atommask(i) == 0) cycle - zi = ((cn(i)**cnexp)*float(at(i))) - do j = 1,i-1 - if (atommask(j) == 0) cycle - r = (xyz(1,i)-xyz(1,j))**2 & - & +(xyz(2,i)-xyz(2,j))**2 & - & +(xyz(3,i)-xyz(3,j))**2 - r = sqrt(r) - zj = ((cn(j)**cnexp)*float(at(j))) - dum = (zi*zj) - dum = dum/(r**rexp) - !dum = dum * exp( -0.5_wp * (r-5)**2) - ecoulomb = ecoulomb+dum - end do - end do - !>-- normalization to the number of included(!) atoms - dum = float(sum(atommask)) - ecoulomb = ecoulomb/(dum**natexp) - ecoulomb = ecoulomb*627.5095 - return -end subroutine calc_ecoul - -!=========================================================================================! - -subroutine cregen_CRE_2(ch,env,nat,nall,at,xyz,comments,nallout,group) -!************************************************************** -!* subroutine cregen_CRE_2 -!* an alternaitve version to the above routine to be used -!* with an unsorted (small) ensemble to determine -!* duplicate groups, but nothing is sorted! -!* On Input: ch - printout channel -!* nat - number of atoms -!* nall - number of structure in ensemble -!* at - atom types -!* xyz - Cartesian coordinates -!* comments - commentary lines containing the energy -!* On Output: nallout - new total number of structures -!* group - to which group every structure belongs -!************************************************************** - use crest_parameters,id => dp - use crest_data - use strucrd - use ls_rmsd - use axis_module - use utilities - implicit none - type(systemdata) :: env - integer,intent(in) :: ch - integer,intent(in) :: nat - integer,intent(in) :: nall - integer,intent(in) :: at(nat) - real(wp),intent(inout) :: xyz(3,nat,nall) - integer,intent(inout) :: group(0:nall) - character(len=*) :: comments(nall) - integer :: nallout - logical :: enantio = .true. !> check for enantiomers? - - !>--- float data - real(wp) :: ewin,rthr,bthr,pthr,ethr,athr,T - !>--- energy data - real(wp),allocatable :: er(:) - integer,allocatable :: orderref(:),order(:) - real(wp) :: de - !>--- dummy structure data - real(wp),allocatable :: c0(:,:),c1(:,:) - real(wp),allocatable :: c0h(:,:),c1h(:,:) - integer,allocatable :: maskheavy(:) - integer,allocatable :: at0(:) - integer :: nat0 - real(wp),allocatable :: rot(:,:) - real(wp) :: bdum - !>--- RMSD data - real(sp),allocatable :: rmat(:) !> SINGLE PRECISION - real(wp) :: rdum,dr,rdum2 - real(wp),allocatable :: gdum(:,:),Udum(:,:),xdum(:),ydum(:) !> dummy tensors - integer(id) :: klong - integer(id),allocatable :: rmap1(:) - integer,allocatable :: rmap2(:) - logical :: l1,l2,l3 - !>--- CRE comparison data - integer,allocatable :: double(:) - logical :: equalrotaniso !> this is a function - logical,allocatable :: mask(:) - real(wp),allocatable :: enuc(:) - real(wp) :: couthr - real(wp),allocatable :: ecoul(:) - real(wp) :: r - integer :: i,j,k,natnoh - logical :: heavy - -!>--- set parameters - call cregen_filldata1(env,ewin,rthr,ethr,bthr,athr,pthr,T,couthr) - if (env%entropic) enantio = .false. - heavy = env%heavyrmsd - allocate (rot(3,nall)) - -!>--- get energies from the comment line - allocate (er(nall)) - allocate (orderref(nall),order(nall)) - do i = 1,nall - er(i) = grepenergy(comments(i)) - orderref(i) = i - end do - -!>--- get dummy structure memory space - allocate (c0(3,nat),c1(3,nat),at0(nat)) - at0 = at - nat0 = nat - -!>--- transform the coordinates to CMA and get rot.constants - do i = 1,nall - call axis(nat,at,xyz(:,:,i)) !>-- all coordinates to CMA - c1(:,:) = xyz(:,:,i) - call axis(nat0,at0,c1,rot(1:3,i),bdum) !>-- B0 in MHz - end do - -!>--- RMSD part - allocate (double(nall),source=0) - !========================================================! - !>-- crucial point: rmat is huge. VERY huge, potentially. - ! use only for small ensembles! - klong = nall - klong = klong*(nall+1) - klong = klong/2 - !write(*,*) nall,klong - allocate (rmat(klong),source=0.0_sp) - allocate (gdum(3,3),Udum(3,3),xdum(3),ydum(3)) - !>-- begin calculation of RMSDs - klong = 0 - write (stdout,'(a)',advance='no') 'CREGEN> running RMSDs ...' - flush (stdout) - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - if (.not.heavy) then !really, the regular case - do i = 1,nall - c0(1:3,1:nat) = xyz(1:3,1:nat,i) -!$OMP PARALLEL PRIVATE ( j,klong,c1,xdum,ydum,Udum,gdum,rdum,rdum2,de) & -!$OMP SHARED ( i,c0,rmat,nat,xyz,rmap1,rmap2,er,ethr,enantio) -!$OMP DO - do j = 1,i-1 - de = (er(i)-er(j))*autokcal - if (de .lt. ethr) then - c1(1:3,1:nat) = xyz(1:3,1:nat,j) - call rmsd(nat,c0,c1,0,Udum,xdum,ydum,rdum,.false.,gdum) ! all atoms - if (enantio) then !also check for enantiomer by inverting a coordinate - c1(1,:) = -c1(1,:) - call rmsd(nat,c0,c1,0,Udum,xdum,ydum,rdum2,.false.,gdum) ! all atoms - else - rdum2 = rdum - end if - !klong=linr(rmap1(i),rmap2(i),j) - klong = lina(i,j) - rmat(klong) = real(min(rdum,rdum2),sp) - end if - end do -!$OMP END DO -!$OMP END PARALLEL - end do - else !heavy atom case - natnoh = nat-counth(nat,at) - allocate (c0h(3,natnoh),c1h(3,natnoh),source=0.0_wp) - allocate (maskheavy(nat),source=0) - call heavymask(nat,at,maskheavy) - write (*,*) 'doing heavy atom rmsds with ',natnoh,' atoms' - do i = 1,nall - call maskedxyz2(nat,natnoh,xyz(:,:,i),c0h,maskheavy) -!$OMP PARALLEL PRIVATE ( j,klong,c1h,xdum,ydum,Udum,gdum,rdum,rdum2,de) & -!$OMP SHARED ( i,c0h,rmat,nat,xyz,rmap1,rmap2,er,ethr,enantio,maskheavy) -!$OMP DO - do j = 1,i-1 - de = (er(i)-er(j))*autokcal - if (de .lt. ethr) then - call maskedxyz2(nat,natnoh,xyz(:,:,j),c1h,maskheavy) - call rmsd(natnoh,c0h,c1h,0,Udum,xdum,ydum,rdum,.false.,gdum) ! all atoms - if (enantio) then !also check for enantiomer by inverting a coordinate - c1h(1,:) = -c1h(1,:) - call rmsd(natnoh,c0h,c1h,0,Udum,xdum,ydum,rdum2,.false.,gdum) ! all atoms - else - rdum2 = rdum - end if - !klong=linr(rmap1(i),rmap2(i),j) - klong = lina(i,j) - rmat(klong) = real(min(rdum,rdum2),sp) - end if - end do -!$OMP END DO -!$OMP END PARALLEL - end do - deallocate (maskheavy,c1h,c0h) + if (nosort) then + write (ch,'(1x,a,i10)') 'number of removed rotamers :', (nallout-k) + nallout = k + write (ch,'(1x,a,i10)') 'total number unique points remaining :',nallout end if - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - write (stdout,'(1x,a)') 'done.' - !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! -!>-- Now, with the RMSDs and rotational constants we can kick out duplicates - do i = 1,nall - do j = 1,i-1 - !>-- only check for structures in energy range - de = (er(i)-er(j))*autokcal - if (de .lt. ethr) then - !klong= linr(rmap1(i),rmap2(i),j) - klong = lina(i,j) - dr = rmat(klong) - else - cycle - end if - !>-- very small RMSD --> same structure - if (dr .lt. rthr) then - double(i) = j !>-- "i" is the same structure as "j" - !>-- slightly larger RMSD, but same rot. constants --> same structure - elseif (dr .lt. 2.0_wp*rthr) then - l1 = equalrotaniso(i,j,nall,rot,0.5d0*bthr,env%bthrmax,env%bthrshift) - if (l1) then - double(i) = j !>-- "i" is the same structure as "j" - end if - end if - end do - !>-- find the original reference. k is a dummy variable - call backtrack(double,i,k) - end do - deallocate (c0,c1,at0) - deallocate (ydum,xdum,Udum,gdum) - deallocate (rmat) - -!>-- for ENSO write a file with duplicate info (if required) - call enso_duplicates(env,nall,double) - -!>-- count how many duplicates we have found - allocate (mask(nall)) - mask(:) = double(:) .ne. 0 - k = count(mask,1) - nallout = nall-k - deallocate (mask) - write (ch,'(1x,a,i10)') 'number of doubles removed by rot/RMSD:',k - !write(ch,*)'total number unique points remaining :',nallout - -!>-- finally, determine conformer groups and their rotamers - allocate (c1(3,nat)) - allocate (enuc(nall)) - do k = 1,nall - c1(1:3,1:nat) = xyz(1:3,1:nat,k) - enuc(k) = 0.0_wp - do i = 1,nat-1 - do j = i+1,nat - r = (c1(1,i)-c1(1,j))**2 & - & +(c1(2,i)-c1(2,j))**2 & - & +(c1(3,i)-c1(3,j))**2+1.d-12 - enuc(k) = enuc(k)+at(i)*at(j)/r - end do - end do - end do -!>-- check energy, rot. const. and nuclear permutation - SORTI: do i = 1,nall - SORTJ: do j = 1,i-1 - !>-- energy difference - de = (er(i)-er(j))*autokcal - l3 = double(j) .eq. 0 - if (.not.l3) cycle - if (abs(de) .lt. ethr) then - !>-- rotational constant difference - l1 = equalrotaniso(i,j,nall,rot,bthr,env%bthrmax,env%bthrshift) - !>-- nuclear permutation - l2 = 2.0d0*abs(enuc(i)-enuc(j))/(enuc(i)+enuc(j)) .lt. 1.d-3 - if (l1.and.l2.and.l3) then - double(i) = j !>-- "i" is a rotamer of "j" - call backtrack(double,i,k) - cycle SORTI - end if - end if - end do SORTJ - end do SORTI -!>-- assign conformer groups - k = 0 - group = 0 - do i = 1,nall - if (double(i) .eq. 0) then - k = k+1 - group(i) = k - else - j = double(i) - group(i) = group(j) - end if - end do - group(0) = k !>-- total number of groups - write (ch,'(1x,a,i10)') 'number of removed rotamers :', (nallout-k) - nallout = k - write (ch,'(1x,a,i10)') 'total number unique points remaining :',nallout deallocate (enuc,c1,double) deallocate (order,orderref) @@ -1824,7 +1499,7 @@ subroutine heavymask(nat,at,mask) end do return end subroutine heavymask -end subroutine cregen_CRE_2 +end subroutine cregen_CRE !=========================================================================================! @@ -2926,7 +2601,7 @@ subroutine cregen_pr3(ch,infile,nall,comments) 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,a,a)') ' written to file <',trim(infile),'>' write (ch,*) write (ch,'(a10,4x,a15,a25)') 'structure','ΔE(kcal/mol)','Etot(Eh)' !write (ch,'('' structure ΔE(kcal/mol) Etot(Eh)'')') diff --git a/src/crest_main.f90 b/src/crest_main.f90 index b5446ce7..f8da3347 100644 --- a/src/crest_main.f90 +++ b/src/crest_main.f90 @@ -350,6 +350,10 @@ program CREST call cs_shutdown(io) end block +!=========================================================================================! +!> one final cleanup + call custom_cleanup(env) + !=========================================================================================! !> Evaluate and print timings call eval_timer(tim) diff --git a/src/strucreader.f90 b/src/strucreader.f90 index b2272e71..12ef7534 100644 --- a/src/strucreader.f90 +++ b/src/strucreader.f90 @@ -127,6 +127,9 @@ module strucrd module procedure wrensemble_conf module procedure wrensemble_conf_energy module procedure wrensemble_conf_energy_comment + + module procedure wrensemble_coord_name + module procedure wrensemble_coord_channel end interface wrensemble public :: pdbdata @@ -496,7 +499,7 @@ end subroutine ensemble_strucskip ! can have a diferent number and order of atoms. ! version 2 does not read energies !=================================================================! - subroutine rdensemble_mixed2(fname,natmax,nall,nats,ats,xyz) + subroutine rdensemble_mixed2(fname,natmax,nall,nats,ats,xyz,comments) implicit none character(len=*),intent(in) :: fname integer,intent(in) :: natmax @@ -504,6 +507,7 @@ subroutine rdensemble_mixed2(fname,natmax,nall,nats,ats,xyz) integer :: nats(nall) integer :: ats(natmax,nall) real(wp) :: xyz(3,natmax,nall) + character(len=*) :: comments(nall) integer :: i,j,k,ich,io logical :: ex integer :: dum @@ -517,6 +521,7 @@ subroutine rdensemble_mixed2(fname,natmax,nall,nats,ats,xyz) nats(i) = dum read (ich,'(a)',iostat=io) line if (io < 0) exit + comments(i) = trim(line) do j = 1,dum read (ich,'(a)',iostat=io) line if (io < 0) exit @@ -552,35 +557,16 @@ subroutine rdensemble_coord_type(fname,nall,structures) integer,allocatable :: at(:) integer,allocatable :: ats(:,:) real(wp),allocatable :: eread(:) + character(len=512),allocatable :: comments(:) integer :: i,j,k,ich,io,nat_i logical :: ex,multiple_sizes call rdensembleparam(fname,nat,nall,multiple_sizes) -! if(.not.multiple_sizes)then -! !>--- 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) + allocate(comments(nall)) + call rdensemble_mixed2(fname,nat,nall,nats,ats,xyz,comments) !>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<--- Important: coord types must be in Bohrs xyz = xyz/bohr @@ -592,11 +578,12 @@ subroutine rdensemble_coord_type(fname,nall,structures) structures(i)%at(:) = ats(1:nat_i,i) allocate (structures(i)%xyz(3,nat_i)) structures(i)%xyz(:,:) = xyz(1:3,1:nat_i,i) + eread(i) = grepenergy(comments(i)) structures(i)%energy = eread(i) + structures(i)%comment = trim(comments(i)) end do - deallocate(eread,nats,ats,xyz) - ! endif + deallocate(comments,eread,nats,ats,xyz) end subroutine rdensemble_coord_type !=================================================================! @@ -685,6 +672,35 @@ subroutine write_ensemble(self,fname) return end subroutine write_ensemble + + subroutine wrensemble_coord_name(fname,nall,structures) + implicit none + character(len=*),intent(in) :: fname + integer,intent(in) :: nall + type(coord) :: structures(nall) + integer :: ich,i + open (newunit=ich,file=fname,status='replace') + do i = 1,nall + call structures(i)%append(ich) + end do + close (ich) + return + end subroutine wrensemble_coord_name + + + subroutine wrensemble_coord_channel(ich,nall,structures) + implicit none + integer,intent(in) :: ich + integer,intent(in) :: nall + type(coord) :: structures(nall) + integer :: i + do i = 1,nall + call structures(i)%append(ich) + end do + return + end subroutine wrensemble_coord_channel + + !==================================================================! ! subroutine deallocate_ensembletype ! is used to clear memory for the ensemble type @@ -2193,21 +2209,33 @@ function grepenergy(line) character(len=*),intent(in) :: line real(wp) :: energy character(len=:),allocatable :: atmp - integer :: i,io + integer :: i,io,k atmp = trim(line) energy = 0.0_wp - !> assumes that the first float in the line is the energy - do i = 1,len_trim(atmp) - if (len_trim(atmp) .lt. 1) exit + if(index(atmp,'energy=').ne.0)then + k=index(atmp,'energy=') + atmp=atmp(k+7:) read (atmp,*,iostat=io) energy - if (io > 0) then - atmp = atmp(2:) - atmp = adjustl(atmp) - cycle - else - exit - end if - end do + if(io.ne.0) energy=0.0_wp + else if(index(atmp,'energy:').ne.0)then + k=index(atmp,'energy:') + atmp=atmp(k+7:) + read (atmp,*,iostat=io) energy + if(io.ne.0) energy=0.0_wp + else + !> assumes that the first float in the line is the energy + do i = 1,len_trim(atmp) + if (len_trim(atmp) .lt. 1) exit + read (atmp,*,iostat=io) energy + if (io > 0) then + atmp = atmp(2:) + atmp = adjustl(atmp) + cycle + else + exit + end if + end do + endif grepenergy = energy return end function grepenergy diff --git a/subprojects/gfnff b/subprojects/gfnff index a49f529d..42963235 160000 --- a/subprojects/gfnff +++ b/subprojects/gfnff @@ -1 +1 @@ -Subproject commit a49f529d3b915350e36ae7743e761b6ac48b817a +Subproject commit 42963235cba66f81575f17b2dba8be7acf2a440d diff --git a/subprojects/tblite b/subprojects/tblite index 82c3dcb9..1665c595 160000 --- a/subprojects/tblite +++ b/subprojects/tblite @@ -1 +1 @@ -Subproject commit 82c3dcb9f76f1b45639b7a17c27e7ba40fffe472 +Subproject commit 1665c5959a588a3b2388b420baca6e9198773cbd