Maintenance (crest-lab#63)

Some additions for maintenance: - optional routine for additional OH orientation sampling (-hflip/-noflip) - improved memory handling for topology detection - added routine for reading a file with atomic charges for GFN-FF calculations (-charges ) Signed-off-by: Philipp Pracht <[email protected]>
pprcht · Jun 29, 2021 · 726ead1 · 726ead1
1 parent a468a58
commit 726ead1
Show file tree

Hide file tree

Showing 23 changed files with 483 additions and 116 deletions.
diff --git a/.gitignore b/.gitignore
@@ -1,14 +1,8 @@
 *.o
 *.mod
-backups/
-old_source/
 *.tgz
-version_pack.sh
 github_bin/
-snippits
-stefan-source/
-schick/
 build_majestix
-crest-build-files.tar.xz
 build_commands
 _build*
+src/crest
diff --git a/src/acidbase.f90 b/src/acidbase.f90
@@ -84,7 +84,7 @@ subroutine acidbase(env,acidfile,basefile,acidchrg,verbose,keepdir,dE, &
       eatb=0.0_wp
       gsa=0.0_wp
       grrhoa=0.0_wp
-      call simpletopo_mol(mola,acid,.false.)
+      call simpletopo_mol(mola,acid,.false.,.false.)
       r = makedir(adir)
       call chdir(adir)
       call mola%write('coord')
@@ -104,7 +104,7 @@ subroutine acidbase(env,acidfile,basefile,acidchrg,verbose,keepdir,dE, &
       ebtb=0.0_wp
       gsb=0.0_wp
       grrhob=0.0_wp
-      call simpletopo_mol(molb,base,.false.)
+      call simpletopo_mol(molb,base,.false.,.false.)
       r = makedir(bdir)
       call chdir(bdir)
       call molb%write('coord')

diff --git a/src/bondconstraint.f90 b/src/bondconstraint.f90
@@ -33,7 +33,7 @@ subroutine autoBondConstraint(filename,forceconstant,wbofile)
      real(wp) :: forceconstant
      type(zmolecule) :: zmol    
    !--- get topology
-     call simpletopo_file(filename,zmol,.false.,wbofile)
+     call simpletopo_file(filename,zmol,.false.,.false.,wbofile)
    !--- get bond matrix and write "bondlengths" file
      call getbmat(zmol,1,forceconstant)
      call zmol%deallocate()
@@ -49,7 +49,7 @@ subroutine autoBondConstraint_withEZ(filename,forceconstant,wbofile)
      real(wp) :: forceconstant
      type(zmolecule) :: zmol    
    !--- get topology
-     call simpletopo_file(filename,zmol,.false.,wbofile)
+     call simpletopo_file(filename,zmol,.false.,.false.,wbofile)
    !--- get bond matrix and write "bondlengths" file
      call getbmat(zmol,5,forceconstant)
      call zmol%deallocate()
@@ -73,7 +73,7 @@ subroutine autoMetalConstraint(filename,forceconstant,wbofile)
      real(wp) :: forceconstant
      type(zmolecule) :: zmol
    !--- get topology
-     call simpletopo_file(filename,zmol,.false.,wbofile)
+     call simpletopo_file(filename,zmol,.false.,.false.,wbofile)
    !--- get bond matrix and write "bondlengths" file
      call getbmat(zmol,2,forceconstant)
      call zmol%deallocate()
@@ -97,7 +97,7 @@ subroutine autoHeavyConstraint(filename,forceconstant,wbofile)
      real(wp) :: forceconstant
      type(zmolecule) :: zmol
    !--- get topology
-     call simpletopo_file(filename,zmol,.false.,'')
+     call simpletopo_file(filename,zmol,.false.,.false.,'')
    !--- get bond matrix and write "bondlengths" file
      call getbmat(zmol,3,forceconstant)
      call zmol%deallocate()
@@ -122,7 +122,7 @@ subroutine autoHydrogenConstraint(filename,forceconstant,wbofile)
      real(wp) :: forceconstant
      type(zmolecule) :: zmol
    !--- get topology
-     call simpletopo_file(filename,zmol,.false.,'')
+     call simpletopo_file(filename,zmol,.false.,.false.,'')
    !--- get bond matrix and write "bondlengths" file
      call getbmat(zmol,4,forceconstant)
      call zmol%deallocate()

diff --git a/src/ccegen.f90 b/src/ccegen.f90
@@ -153,14 +153,13 @@ subroutine CCEGEN(env,pr,fname)
        fraclimit = 0.25d0 !if autolimit=true, 1/4 of the ensemble is taken
 
    !--- 1. topology for reference strucuture
-      !call simpletopo_file('coord',zmol,pr)
       if(env%wbotopo)then
          env%wbofile='wbo'
       else   
          env%wbofile='none given' 
       endif
       zens%xyz = zens%xyz/bohr  !ANG to Bohr for topo
-      call simpletopo(zens%nat,zens%at,zens%xyz,zmol,pr,env%wbofile)
+      call simpletopo(zens%nat,zens%at,zens%xyz,zmol,pr,.false.,env%wbofile)
       zens%xyz = zens%xyz*bohr  !Bohr to ANG
       allocate(inc(zmol%nat), source = 0)
 

diff --git a/src/classes.f90 b/src/classes.f90
@@ -243,6 +243,9 @@ module crest_data
        real(wp),allocatable :: xyz(:,:)
        integer :: ntopo
        integer,allocatable :: topo(:)
+       real(wp),allocatable :: charges(:)
+    contains
+      procedure :: rdcharges => read_charges
    end type refdata
 
 !-----------------------------------------------------------------------------------------------------
@@ -321,7 +324,8 @@ module crest_data
       character(len=:),allocatable :: inputcoords
       character(len=:),allocatable :: wbofile
       character(len=:),allocatable :: atlist
-
+      character(len=:),allocatable :: chargesfilename
+
      !--- METADYN data
       real(wp) :: hmass
       real(wp) :: mdtemp
@@ -386,6 +390,10 @@ module crest_data
       character(len=:),allocatable :: pcmeasure
       integer :: clustlev = 0    ! clustering level
 
+    !--- additional structure generation settings
+      logical :: doOHflip = .false.
+      integer :: maxflip = 1000  
+
     !--- external RMSD bias to optimizations
       character(len=:),allocatable :: biasfile
       real(wp) :: rthr2 = 0.3_wp    ! Discard all structures with a bias smaller than this
@@ -419,6 +427,7 @@ module crest_data
       logical :: cluster = .false. ! perform a clustering analysis        
       logical :: checktopo = .true.  !perform topolgy check in CREGEN
       logical :: checkiso  = .false. !perform E/Z isomerization check in CREGEN
+      logical :: chargesfile = .false. !use a given charges file for gfnff
       logical :: compareens        ! try to correlate 2 given Ensemble files
       logical :: confgo            ! perform only the CREGEN routine ?
       logical :: doNMR             ! determine NMR equivalencies in CREGEN ?
@@ -464,6 +473,7 @@ module crest_data
       logical :: relax =.false.    ! was the --relax function used for protonation site search?
       logical :: restartopt        ! restart in the second step of the multilevel opt (V2) ?
       logical :: reweight=.false.  ! reweight structures on the fly after optimizations (i.e. do SPs)?
+      logical :: riso=.false.      ! take only isomers in reactor mode
       logical :: rotamermds        ! do additional MDs after second  multilevel OPT step in V2 ?
       logical :: scallen           ! scale the automatically determined MD length by some factor?
       logical :: scratch           ! use scratch directory
@@ -705,7 +715,7 @@ subroutine wrtCHRG(self,dir)
     class(systemdata) :: self
     character(len=*) :: dir    
     character(len=:),allocatable :: path
-    integer :: ich,k
+    integer :: ich,k,i
     k = len_trim(dir)
     if(self%chrg.ne.0)then
        if(k>0)then 
@@ -727,8 +737,49 @@ subroutine wrtCHRG(self,dir)
        write(ich,*) self%uhf
        close(ich)
     endif  
+    if(self%chargesfile .and. allocated(self%ref%charges))then
+       if(k>0)then 
+       path=trim(dir)//'/'//'charges'
+       else
+       path='charges'
+       endif
+       open(newunit=ich,file=path)
+       do i=1,self%ref%nat
+       write(ich,'(1x,f16.8)') self%ref%charges(i)
+       enddo
+       close(ich)
+    endif
+    return
+end subroutine wrtCHRG 
+!------------------------------------------------------------------------------------------------------
+! read atomic charges from a file (one line per atom)
+subroutine read_charges(self,chargesfilename,totchrg)
+    implicit none
+    class(refdata) :: self
+    character(len=*) :: chargesfilename
+    integer :: ich,io,i
+    real(wp) :: dum,tot
+    integer :: totchrg
+    if(allocated(self%charges))deallocate(self%charges)
+    if(self%nat>0)then
+        allocate(self%charges(self%nat),source=0.0_wp)
+        open(newunit=ich,file=chargesfilename)
+        do i=1,self%nat
+         read(ich,*,iostat=io) dum
+         if(io==0)then
+           self%charges(i) = dum
+         endif
+        enddo
+        close(io)
+    endif
+    tot=0.0_wp
+    do i=1,self%nat
+      tot=tot+self%charges(i)
+    enddo
+    totchrg = nint(tot)
     return
-end subroutine wrtCHRG    
+end subroutine read_charges
+
 
 !------------------------------------------------------------------------------------------------------
 function optlevflag(optlev) result(flag)

diff --git a/src/confparse.f90 b/src/confparse.f90
@@ -41,7 +41,7 @@ subroutine parseflags(env,arg,nra)  !FOR THE CONFSCRIPT STANDALONE
 
       character(len=512) :: atmp,btmp
       character(len=:),allocatable :: ctmp,dtmp
-      integer :: i,j,k,l,io,ich
+      integer :: i,j,k,l,io,ich,idum
       integer*8 :: w,v
       real(wp) :: rdum
 
@@ -340,6 +340,7 @@ subroutine parseflags(env,arg,nra)  !FOR THE CONFSCRIPT STANDALONE
           case( '-mdopt','-purge' )                                           ! MDOPT
                env%crestver = crest_mdopt
                atmp=''
+               env%preopt=.false.
                env%ensemblename='none selected'
                if(nra.ge.(i+1))atmp=adjustl(arg(i+1))
                if((atmp(1:1)/='-').and.(len_trim(atmp).ge.1))then
@@ -607,7 +608,7 @@ subroutine parseflags(env,arg,nra)  !FOR THE CONFSCRIPT STANDALONE
                 call redo_extrapol(env,ctmp,0)    
                 endif
                 stop
-           case('-SANDBOX' )
+              case('-SANDBOX' )
             !--- IMPLEMENT HERE WHATEVER YOU LIKE, FOR TESTING
 
             !-----
@@ -895,6 +896,8 @@ subroutine parseflags(env,arg,nra)  !FOR THE CONFSCRIPT STANDALONE
                    env%preactormtd=.true.
                 case( '-fragopt' ) 
                    env%restartopt=.true.
+                case( '-iso' )
+                   env%riso=.true.   
                 case default
                     continue
                end select RCTR
@@ -985,6 +988,22 @@ subroutine parseflags(env,arg,nra)  !FOR THE CONFSCRIPT STANDALONE
                  else
                  write(*,'(2x,a,a)') argument,' : energy reweighting'
                  endif    
+             case('-charges') !read charges from file for GFN-FF calcs.
+                ctmp=trim(arg(i+1))
+                if((len_trim(ctmp)<1).or.(ctmp(1:1)=='-'))then
+                    ctmp='charges'
+                endif
+                inquire(file=ctmp,exist=ex)
+                if(ex)then
+                    env%chargesfilename=ctmp
+                    env%chargesfile=.true.
+                    write(*,'(2x,a,a,a)') '-charges: file <',trim(ctmp),'> used for atomic charges'
+                    call env%ref%rdcharges(env%chargesfilename,idum)
+                    if(idum.ne.env%chrg)then
+                    write(*,'(12x,a,i0)') 'with total summed up molecular charge: ',idum
+                    env%chrg = idum
+                    endif
+                endif
              case( '-dscal','-dispscal','-dscal_global','-dispscal_global' )
                  env%cts%dispscal_md = .true.
                  if(index(argument,'_global').ne.0)then
@@ -1183,6 +1202,15 @@ subroutine parseflags(env,arg,nra)  !FOR THE CONFSCRIPT STANDALONE
                 read(arg(i+1),*,iostat=io) rdum
                 if(io==0) env%kshift = rdum   
                 env%kshiftnum = 1
+             case( '-hflip' )
+                env%doOHflip = .true.
+             case( '-noflip' )
+                env%doOHflip = .false.
+             case( '-maxflip' )
+                read(arg(i+1),*,iostat=io) rdum
+                if(io==0 .and. (index(arg(i+1),'-').eq.0))then
+                  env%maxflip = nint(rdum)
+                endif
 !======================================================================================================!
 !------ flags for parallelization / disk space
 !======================================================================================================!
@@ -1797,7 +1825,7 @@ subroutine parseRC2(env,bondconst)
       endif
 
       if(ex1)then
-         write(*,'(1x,a,a,a)') '"',trim(env%constraints),'" file present.'
+         write(*,'(/,1x,a,a,a)') '<',trim(env%constraints),'> file present.'
          cts%used=.true.
       else
          cts%used=.false.
@@ -1981,7 +2009,7 @@ subroutine inputcoords(env,arg)
 
     !--- for protonation/deprotonation applications get ref. number of fragments
     !--- also get some other structure based info
-    call simpletopo_file('coord',zmol,.false.,'')
+    call simpletopo_file('coord',zmol,.false.,.false.,'')
     env%ptb%nfrag = zmol%nfrag
     call zmol%deallocate()
 

diff --git a/src/confscript2_main.f90 b/src/confscript2_main.f90
@@ -202,6 +202,24 @@ subroutine confscript2i(env,tim)
        call remaining_in(atmp,ewin,nallout) !--- remaining number of structures
        write(*,*)
 !====================================================================!
+!---- (Optional) sampling of additional XH positions
+      if(env%doOHflip)then
+          call XHorient(env,conformerfile)
+          inquire(file='oh_ensemble.xyz',exist=ex)
+          if(ex)then
+            call checkname_xyz(crefile,atmp,btmp)
+            call appendto('oh_ensemble.xyz',atmp)       
+            call remove('oh_ensemble.xyz')
+            if(.not.env%entropic .and. env%crestver.ne.22)then
+              call newcregen(env,0)
+            else
+              call newcregen(env,2)
+            endif
+            call remaining_in(btmp,ewin,nallout)
+            write(*,*)
+          endif
+      endif
+!====================================================================!
 !---- Perform additional MDs on the lowest conformers
       if(env%rotamermds)then
           call tim%start(4,'MD ')