1! This file is part of xtb. 2! 3! Copyright (C) 2019-2020 Stefan Grimme 4! 5! xtb is free software: you can redistribute it and/or modify it under 6! the terms of the GNU Lesser General Public License as published by 7! the Free Software Foundation, either version 3 of the License, or 8! (at your option) any later version. 9! 10! xtb is distributed in the hope that it will be useful, 11! but WITHOUT ANY WARRANTY; without even the implied warranty of 12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13! GNU Lesser General Public License for more details. 14! 15! You should have received a copy of the GNU Lesser General Public License 16! along with xtb. If not, see <https://www.gnu.org/licenses/>. 17 18module xtb_gfnff_setup 19 use xtb_gfnff_ini, only : gfnff_ini 20 use xtb_gfnff_data, only : TGFFData 21 use xtb_gfnff_topology, only : TGFFTopology 22 use xtb_gfnff_generator, only : TGFFGenerator 23 implicit none 24 private 25 public :: gfnff_setup, gfnff_input 26 27contains 28 29subroutine gfnff_setup(env,verbose,restart,mol,gen,param,topo,accuracy,version) 30 use xtb_restart 31 use xtb_type_environment, only : TEnvironment 32 use xtb_type_molecule, only : TMolecule 33 use xtb_gfnff_param, only : ini, gfnff_set_param 34 use xtb_setparam, only : ichrg, dispscale 35 implicit none 36 character(len=*), parameter :: source = 'gfnff_setup' 37! Dummy 38 !integer,intent(in) :: ich 39 type(TGFFTopology), intent(inout) :: topo 40 type(TGFFGenerator), intent(inout) :: gen 41 type(TGFFData), intent(inout) :: param 42 integer,intent(in) :: version 43 logical,intent(in) :: restart 44 logical,intent(in) :: verbose 45 real(wp),intent(in) :: accuracy 46 type(TMolecule) :: mol 47 type(TEnvironment), intent(inout) :: env 48! Stack 49 logical :: ex 50 logical :: success 51 logical :: exitRun 52 53 call gfnff_input(env, mol, topo) 54 call env%check(exitRun) 55 if (exitRun) then 56 call env%error("Failed to prepare topology from geometry input", source) 57 return 58 end if 59 60 call gfnff_set_param(mol%n, gen, param) 61 param%dispscale = dispscale 62 if (restart) then 63 inquire(file='gfnff_topo', exist=ex) 64 if (ex) then 65 call read_restart_gff(env,'gfnff_topo',mol%n,version,success,.true.,topo) 66 if (success) then 67 write(env%unit,'(10x,"GFN-FF topology read from file successfully!")') 68 return 69 else 70 call env%warning("Could not read topology file.", source) 71 call env%check(exitRun) 72 if (exitRun) then 73 return 74 end if 75 76 end if 77 end if 78 end if 79 80 call gfnff_ini(env,verbose,ini,mol,ichrg,gen,param,topo,accuracy) 81 82 call env%check(exitRun) 83 if (exitRun) then 84 call env%error("Failed to generate topology", source) 85 return 86 end if 87 88 if (.not.mol%struc%two_dimensional) then 89 call write_restart_gff(env,'gfnff_topo',mol%n,version,topo) 90 end if 91 92end subroutine gfnff_setup 93 94subroutine gfnff_input(env, mol, topo) 95 use xtb_mctc_accuracy, only : wp 96 use xtb_type_environment, only : TEnvironment 97 use xtb_type_molecule 98 use xtb_mctc_filetypes, only : fileType 99 use xtb_gfnff_param 100 use xtb_setparam, only : ichrg 101 implicit none 102 ! Dummy 103 type(TMolecule),intent(in) :: mol 104 type(TEnvironment), intent(inout) :: env 105 type(TGFFTopology), intent(inout) :: topo 106 ! Stack 107 integer :: i,j,k 108 integer :: ni 109 integer :: ns 110 integer :: nf 111 integer :: ich 112 integer :: iatom 113 integer :: iresidue 114 integer :: ifrag 115 integer :: ibond 116 integer :: bond_ij(3) 117 real(wp) :: r 118 real(wp) :: dum1 119 real(wp) :: floats(10) 120 logical :: ex 121 character(len=80) :: atmp 122 character(len=80) :: s(10) 123 integer, allocatable :: rn(:) 124 125 if (.not.allocated(topo%nb)) allocate( topo%nb(20,mol%n), source = 0 ) 126 if (.not.allocated(topo%qfrag)) allocate( topo%qfrag(mol%n), source = 0.0d0 ) 127 if (.not.allocated(topo%fraglist)) allocate( topo%fraglist(mol%n), source = 0 ) 128 129 select case(mol%ftype) 130 !-------------------------------------------------------------------- 131 ! PDB case 132 case(fileType%pdb) 133 ini = .true. 134 ifrag=0 135 allocate(rn(mol%n)) 136 rn(:) = mol%pdb%residue_number 137 do iresidue = minval(rn),maxval(rn) 138 if (any(iresidue .eq. rn)) then 139 ifrag=ifrag+1 140 where(iresidue .eq. rn) topo%fraglist = ifrag 141 end if 142 end do 143 deallocate(rn) 144 topo%nfrag = maxval(topo%fraglist) 145 if (.not.allocated(topo%qpdb)) allocate(topo%qpdb(mol%n)) 146 do iatom=1,mol%n 147 topo%qfrag(topo%fraglist(iatom)) = topo%qfrag(topo%fraglist(iatom)) & 148 & + dble(mol%pdb(iatom)%charge) 149 topo%qpdb(iatom) = mol%pdb(iatom)%charge 150 end do 151 ichrg=idint(sum(topo%qfrag(1:topo%nfrag))) 152 write(env%unit,'(10x,"charge from pdb residues: ",i0)') ichrg 153 !-------------------------------------------------------------------- 154 ! SDF case 155 case(fileType%sdf,fileType%molfile) 156 ini = .false. 157 topo%nb=0 158 topo%nfrag=0 159 do ibond = 1, len(mol%bonds) 160 call mol%bonds%get_item(ibond,bond_ij) 161 i = bond_ij(1) 162 j = bond_ij(2) 163 ni=topo%nb(20,i) 164 ex=.false. 165 do k=1,ni 166 if(topo%nb(k,i).eq.j) then 167 ex=.true. 168 exit 169 endif 170 enddo 171 if(.not.ex)then 172 topo%nb(20,i)=topo%nb(20,i)+1 173 topo%nb(topo%nb(20,i),i)=j 174 topo%nb(20,j)=topo%nb(20,j)+1 175 topo%nb(topo%nb(20,j),j)=i 176 endif 177 end do 178 do i=1,mol%n 179 if(topo%nb(20,i).eq.0)then 180 dum1=1.d+42 181 do j=1,i 182 r=sqrt(sum((mol%xyz(:,i)-mol%xyz(:,j))**2)) 183 if(r.lt.dum1.and.r.gt.0.001)then 184 dum1=r 185 k=j 186 endif 187 enddo 188 topo%nb(20,i)=1 189 topo%nb(1,i)=k 190 endif 191 end do 192 ! initialize qfrag as in the default case 193 topo%qfrag(1)=mol%chrg 194 topo%qfrag(2:mol%n)=0 195 !-------------------------------------------------------------------- 196 ! General case: input = xyz or coord 197 case default 198 if (mol%npbc > 0) then 199 call env%error("Input file format not suitable for GFN-FF!") 200 return 201 end if 202 ini = .true. 203 call open_file(ich,'.CHRG','r') 204 if (ich.ne.-1) then 205 read(ich,'(a)')atmp 206 call close_file(ich) 207 call readline(atmp,floats,s,ns,nf) 208 topo%qfrag(1:nf)=floats(1:nf) 209 ichrg=int(sum(topo%qfrag(1:nf))) 210 topo%qfrag(nf+1:mol%n)=9999 211 else 212 topo%qfrag(1)=mol%chrg 213 topo%qfrag(2:mol%n)=0 214 end if 215 end select 216 217 !------------------------------------------------------------------- 218 219end subroutine gfnff_input 220 221end module xtb_gfnff_setup 222