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