1! This file is part of xtb.
2!
3! Copyright (C) 2019-2020 Sebastian Ehlert
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
18!> Compare and generate topological information from wavefunctions
19module xtb_topology
20   use xtb_mctc_accuracy, only : wp
21   use xtb_mctc_filetypes, only : fileType, generateFileName
22   use xtb_io_writer, only : writeMolecule
23   use xtb_type_molecule, only : TMolecule, len
24   use xtb_type_neighbourlist, only : TNeighbourList, init, resizeNeigh
25   use xtb_type_topology, only : TTopology, len
26   implicit none
27   private
28
29   public :: checkTopology
30   public :: makeBondTopology, compareBondTopology
31
32
33contains
34
35
36subroutine checkTopology(unit, mol, wbo, verbosity)
37
38   integer, intent(in) :: unit
39
40   type(TMolecule), intent(in) :: mol
41
42   real(wp), intent(in) :: wbo(:, :)
43
44   integer, intent(in) :: verbosity
45
46   character(len=:), allocatable :: fname
47   type(TMolecule) :: copy
48   integer :: ich, ftype
49   logical :: match
50
51   ! CT files have some limits:
52   ! bail out if we encounter a structure not suitable for this format
53   if (mol%npbc /= 0 .or. mol%n > 999) return
54
55   copy = mol
56   call makeBondTopology(copy%bonds, mol, wbo)
57
58   call compareBondTopology(unit, mol%bonds, copy%bonds, mol, verbosity, match)
59
60   if (.not.match .and. len(copy%bonds) > 0) then
61      ! prefer molfiles, except the input has been SDF
62      if (mol%ftype /= fileType%sdf) then
63         ftype = fileType%molfile
64      else
65         ftype = mol%ftype
66      end if
67      call generateFileName(fname, 'xtbtopo', '', ftype)
68      call open_file(ich, fname, 'w')
69      if (verbosity > 0 .and. len(mol%bonds) > 0) then
70         write(unit, '("Writing corrected topology to", 1x, a)') fname
71      else
72         write(unit, '("Writing topology from bond orders to", 1x, a)') fname
73      end if
74      call writeMolecule(copy, ich, format=ftype)
75      call close_file(ich)
76   end if
77
78end subroutine checkTopology
79
80
81subroutine makeBondTopology(topo, mol, wbo)
82
83   type(TTopology), intent(out) :: topo
84
85   type(TMolecule), intent(in) :: mol
86
87   real(wp), intent(in) :: wbo(:, :)
88
89   real(wp), parameter :: wbothr = 0.3_wp
90   integer :: iat, jat
91
92   call topo%allocate(size=len(mol), order=3)
93   do iat = 1, len(mol)
94      do jat = 1, iat-1
95         associate(w => wbo(jat, iat))
96            if (any(mol%at(iat) == [6, 7, 8]) .and. any(mol%at(jat) == [6, 7, 8])) then
97               if (w > wbothr .and. w <= 1.2_wp) then
98                  call topo%push_back([jat, iat, 1])
99               else if (w > 1.5_wp .and. w <= 2.3_wp) then
100                  call topo%push_back([jat, iat, 2])
101               else if (w > 2.3_wp) then
102                  call topo%push_back([jat, iat, 3])
103               else if (w > 1.2_wp .and. w <= 1.5_wp) then
104                  call topo%push_back([jat, iat, 4])
105               end if
106            else
107               if (w > wbothr .and. w <= 1.3_wp) then
108                  call topo%push_back([jat, iat, 1])
109               else if (w > 1.3_wp .and. w <= 2.3_wp) then
110                  call topo%push_back([jat, iat, 2])
111               else if (w > 2.3_wp) then
112                  call topo%push_back([jat, iat, 3])
113               end if
114            end if
115         end associate
116      end do
117   end do
118end subroutine makeBondTopology
119
120
121subroutine compareBondTopology(unit, topo1, topo2, mol, verbosity, match)
122
123   integer, intent(in) :: unit
124
125   type(TTopology), intent(in) :: topo1
126
127   type(TTopology), intent(in) :: topo2
128
129   type(TMolecule), intent(in) :: mol
130
131   integer, intent(in) :: verbosity
132
133   logical, intent(out) :: match
134
135   type(TNeighbourList) :: neighList1, neighList2
136   integer :: iat, neigh
137
138   match = len(topo1) == len(topo2)
139
140   if (.not.match .and. verbosity > 0) then
141      write(unit, '("Topologies differ in total number of bonds")')
142   end if
143
144   if (match .or. verbosity > 1) then
145      call topologyToNeighbourList(topo1, neighList1, mol)
146      call topologyToNeighbourList(topo2, neighList2, mol)
147
148      match = match .and. all(neighList1%neighs == neighList2%neighs)
149      if (.not.match .and. verbosity > 0) then
150         write(unit, '("Topologies differ in number of bonds")')
151      end if
152
153      if (match .or. verbosity > 1) then
154         neigh = min(ubound(neighList1%iNeigh, 1), ubound(neighList2%iNeigh, 1))
155         match = match .and. all(neighList1%iNeigh(:neigh, :) &
156            & == neighList2%iNeigh(:neigh, :))
157         if (.not.match .and. verbosity > 0) then
158            write(unit, '("Topologies differ in connectivity")')
159         end if
160
161         if (match .or. verbosity > 1) then
162            match = match .and. all(nint(neighList1%weight(:neigh, :)) &
163               & == nint(neighList2%weight(:neigh, :)))
164            if (.not.match .and. verbosity > 0) then
165               write(unit, '("Topologies differ in bond orders")')
166            end if
167         end if
168      end if
169   end if
170
171end subroutine compareBondTopology
172
173
174subroutine topologyToNeighbourList(topo, neighList, mol)
175
176   type(TTopology), intent(in) :: topo
177
178   type(TMolecule), intent(in) :: mol
179
180   type(TNeighbourList), intent(out) :: neighList
181
182   integer :: iBond(3)
183   integer :: mNeigh
184   integer :: ii, iat, jat
185   real(wp) :: r2
186
187   call init(neighList, len(mol))
188
189   mNeigh = ubound(neighList%iNeigh, dim=1)
190
191   neighList%neighs(:) = 0
192   neighList%iNeigh(:, :) = 0
193   neighList%dist2(:, :) = 0.0_wp
194   neighList%weight(:, :) = 0.0_wp
195   do iat = 1, len(mol)
196      neighList%iNeigh(0, iat) = iat
197      neighList%weight(0, iat) = 1.0_wp
198      neighList%coords(:, iat) = mol%xyz(:, iat)
199      neighList%image(iat) = iat
200      neighList%trans(iat) = 1
201   end do
202
203   do ii = 1, len(topo)
204      call topo%get_item(ii, iBond)
205      iat = min(iBond(1), iBond(2))
206      jat = max(iBond(1), iBond(2))
207      r2 = sum((neighList%coords(:, iat) - neighList%coords(:, jat))**2)
208      neighList%neighs(iat) = neighList%neighs(iat) + 1
209      if (neighList%neighs(iat) > mNeigh) then
210         mNeigh = 2*mNeigh
211         call resizeNeigh(mNeigh, neighList%iNeigh, neighList%dist2, &
212            & neighList%weight)
213      end if
214      neighList%iNeigh(neighList%neighs(iat), iat) = jat
215      neighList%dist2(neighList%neighs(iat), iat) = r2
216      neighList%weight(neighList%neighs(iat), iat) = iBond(3)
217   end do
218
219   call neighList%sort
220
221end subroutine topologyToNeighbourList
222
223
224end module xtb_topology
225