1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  1993  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     ###############################################################
9c     ##                                                           ##
10c     ##  subroutine kopbend  --  out-of-plane bending parameters  ##
11c     ##                                                           ##
12c     ###############################################################
13c
14c
15c     "kopbend" assigns the force constants for out-of-plane bends
16c     at trigonal centers via Wilson-Decius-Cross or Allinger angles;
17c     also processes any new or changed parameter values
18c
19c
20      subroutine kopbend
21      use angbnd
22      use atomid
23      use atoms
24      use couple
25      use fields
26      use inform
27      use iounit
28      use keys
29      use kopbnd
30      use opbend
31      use potent
32      use usage
33      implicit none
34      integer i,j,it
35      integer ia,ib,ic,id
36      integer ita,itb,itc,itd
37      integer nopb,size
38      integer next,number
39      real*8 fopb
40      logical header,done
41      logical, allocatable :: jopb(:)
42      character*4 pa,pb,pc,pd
43      character*4 zero4
44      character*8 zero8
45      character*16 blank,pt
46      character*16 pt0,pt1
47      character*20 keyword
48      character*240 record
49      character*240 string
50c
51c
52c     process keywords containing out-of-plane bend parameters
53c
54      blank = '                '
55      zero4 = '0000'
56      zero8 = '00000000'
57      header = .true.
58      do i = 1, nkey
59         next = 1
60         record = keyline(i)
61         call gettext (record,keyword,next)
62         call upcase (keyword)
63         if (keyword(1:7) .eq. 'OPBEND ') then
64            ia = 0
65            ib = 0
66            ic = 0
67            id = 0
68            fopb = 0.0d0
69            string = record(next:240)
70            read (string,*,err=10,end=10)  ia,ib,ic,id,fopb
71   10       continue
72            size = 4
73            call numeral (ia,pa,size)
74            call numeral (ib,pb,size)
75            call numeral (ic,pc,size)
76            call numeral (id,pd,size)
77            if (ic .le. id) then
78               pt = pa//pb//pc//pd
79            else
80               pt = pa//pb//pd//pc
81            end if
82            if (.not. silent) then
83               if (header) then
84                  header = .false.
85                  write (iout,20)
86   20             format (/,' Additional Out-of-Plane Bend',
87     &                       ' Parameters :',
88     &                    //,5x,'Atom Classes',19x,'K(OPB)',/)
89               end if
90               write (iout,30)  ia,ib,ic,id,fopb
91   30          format (4x,4i4,10x,f12.3)
92            end if
93            size = 4
94            do j = 1, maxnopb
95               if (kopb(j).eq.blank .or. kopb(j).eq.pt) then
96                  kopb(j) = pt
97                  opbn(j) = fopb
98                  goto 50
99               end if
100            end do
101            write (iout,40)
102   40       format (/,' KOPBEND --  Too many Out-of-Plane',
103     &                 ' Angle Bending Parameters')
104            abort = .true.
105   50       continue
106         end if
107      end do
108c
109c     perform dynamic allocation of some global arrays
110c
111      if (allocated(iopb))  deallocate (iopb)
112      if (allocated(opbk))  deallocate (opbk)
113      allocate (iopb(nangle))
114      allocate (opbk(nangle))
115c
116c     use special out-of-plane bend parameter assignment for MMFF
117c
118      if (forcefield .eq. 'MMFF94') then
119         call kopbendm
120         return
121      end if
122c
123c     determine the total number of forcefield parameters
124c
125      nopb = maxnopb
126      do i = maxnopb, 1, -1
127         if (kopb(i) .eq. blank)  nopb = i - 1
128      end do
129c
130c     perform dynamic allocation of some local arrays
131c
132      allocate (jopb(maxclass))
133c
134c     make list of atom classes using out-of-plane bending
135c
136      do i = 1, maxclass
137         jopb(i) = .false.
138      end do
139      do i = 1, maxnopb
140         if (kopb(i) .eq. blank)  goto 60
141         it = number(kopb(i)(5:8))
142         jopb(it) = .true.
143      end do
144   60 continue
145c
146c     assign out-of-plane bending parameters for each angle
147c
148      nopbend = 0
149      if (nopb .ne. 0) then
150         header = .true.
151         do i = 1, nangle
152            ib = iang(2,i)
153            itb = class(ib)
154            if (jopb(itb) .and. n12(ib).eq.3) then
155               ia = iang(1,i)
156               ita = class(ia)
157               ic = iang(3,i)
158               itc = class(ic)
159               id = iang(4,i)
160               itd = class(id)
161               size = 4
162               call numeral (ita,pa,size)
163               call numeral (itb,pb,size)
164               call numeral (itc,pc,size)
165               call numeral (itd,pd,size)
166               if (ita .le. itc) then
167                  pt = pd//pb//pa//pc
168               else
169                  pt = pd//pb//pc//pa
170               end if
171               pt1 = pd//pb//zero8
172               pt0 = zero4//pb//zero8
173               done = .false.
174               do j = 1, nopb
175                  if (kopb(j) .eq. pt) then
176                     nopbend = nopbend + 1
177                     iopb(nopbend) = i
178                     opbk(nopbend) = opbn(j)
179                     done = .true.
180                     goto 70
181                  end if
182               end do
183               do j = 1, nopb
184                  if (kopb(j) .eq. pt1) then
185                     nopbend = nopbend + 1
186                     iopb(nopbend) = i
187                     opbk(nopbend) = opbn(j)
188                     done = .true.
189                     goto 70
190                  end if
191               end do
192               do j = 1, nopb
193                  if (kopb(j) .eq. pt0) then
194                     nopbend = nopbend + 1
195                     iopb(nopbend) = i
196                     opbk(nopbend) = opbn(j)
197                     done = .true.
198                     goto 70
199                  end if
200               end do
201   70          continue
202               if (use_opbend .and. .not.done) then
203                  if (use(ia) .or. use(ib) .or. use(ic) .or. use(id))
204     &               abort = .true.
205                  if (header) then
206                     header = .false.
207                     write (iout,80)
208   80                format (/,' Undefined Out-of-Plane Bend',
209     &                          ' Parameters :',
210     &                       //,' Type',24x,'Atom Names',24x,
211     &                          'Atom Classes',/)
212                  end if
213                  write (iout,90)  id,name(id),ib,name(ib),ia,name(ia),
214     &                             ic,name(ic),itd,itb,ita,itc
215   90             format (' Angle-OP',3x,4(i6,'-',a3),5x,4i5)
216               end if
217            else
218               iang(4,i) = ib
219            end if
220         end do
221      end if
222c
223c     perform deallocation of some local arrays
224c
225      deallocate (jopb)
226c
227c     turn off the out-of-plane bending term if it is not used
228c
229      if (nopbend .eq. 0)  use_opbend = .false.
230      return
231      end
232c
233c
234c     ##################################################################
235c     ##                                                              ##
236c     ##  subroutine kopbendm  --  MMFF out-of-plane bend parameters  ##
237c     ##                                                              ##
238c     ##################################################################
239c
240c
241c     "kopbendm" assigns the force constants for out-of-plane bends
242c     according to the Merck Molecular Force Field (MMFF)
243c
244c
245      subroutine kopbendm
246      use angbnd
247      use atomid
248      use atoms
249      use kopbnd
250      use merck
251      use opbend
252      use potent
253      implicit none
254      integer i,j,m
255      integer nopb,size
256      integer ia,ib,ic,id
257      integer ita,itb,itc,itd
258      integer itta,ittb
259      integer ittc,ittd
260      character*4 pa,pb,pc,pd
261      character*16 blank,pt
262c
263c
264c     determine the total number of forcefield parameters
265c
266      blank = '                '
267      nopb = maxnopb
268      do i = maxnopb, 1, -1
269         if (kopb(i) .eq. blank)  nopb = i - 1
270      end do
271c
272c     assign MMFF out-of-plane bending parameter values
273c
274      nopbend = 0
275      if (nopb .ne. 0) then
276         do i = 1, nangle
277            ia = iang(1,i)
278            ib = iang(2,i)
279            ic = iang(3,i)
280            id = iang(4,i)
281            if (min(ia,ib,ic,id) .gt. 0) then
282               itta = type(ia)
283               ittb = type(ib)
284               ittc = type(ic)
285               ittd = type(id)
286               m = 0
287   10          continue
288               m = m + 1
289               if (m .eq. 1) then
290                  ita = eqclass(itta,1)
291                  itb = eqclass(ittb,1)
292                  itc = eqclass(ittc,1)
293                  itd = eqclass(ittd,1)
294               else if (m .eq. 2) then
295                  ita = eqclass(itta,2)
296                  itb = eqclass(ittb,2)
297                  itc = eqclass(ittc,2)
298                  itd = eqclass(ittd,2)
299               else if (m .eq. 3) then
300                  ita = eqclass(itta,3)
301                  itb = eqclass(ittb,2)
302                  itc = eqclass(ittc,3)
303                  itd = eqclass(ittd,3)
304               else if (m .eq. 4) then
305                  ita = eqclass(itta,4)
306                  itb = eqclass(ittb,2)
307                  itc = eqclass(ittc,4)
308                  itd = eqclass(ittd,4)
309               else if (m .eq. 5) then
310                  ita = eqclass(itta,5)
311                  itb = eqclass(ittb,2)
312                  itc = eqclass(ittc,5)
313                  itd = eqclass(ittd,5)
314               end if
315               if (m .gt. 5) then
316                  nopbend = nopbend + 1
317                  iopb(nopbend) = i
318                  opbk(nopbend) = 0.0d0
319               else
320                  size = 4
321                  call numeral (ita,pa,size)
322                  call numeral (itb,pb,size)
323                  call numeral (itc,pc,size)
324                  call numeral (itd,pd,size)
325                  if (itd.le.ita .and. itd.le.itc) then
326                     if (ita .le. itc) then
327                        pt = pd//pb//pa//pc
328                     else
329                        pt = pd//pb//pc//pa
330                     end if
331                  else if (ita.le.itc .and. ita.le.itd) then
332                     if (itd .le. itc) then
333                        pt = pa//pb//pd//pc
334                     else
335                        pt = pa//pb//pc//pd
336                     end if
337                  else if (itc.le.ita .and. itc.le.itd) then
338                     if (ita .le. itd) then
339                        pt = pc//pb//pa//pd
340                     else
341                        pt = pc//pb//pd//pa
342                     end if
343                  end if
344                  do j = 1, nopb
345                     if (kopb(j) .eq. pt) then
346                        nopbend = nopbend + 1
347                        iopb(nopbend) = i
348                        opbk(nopbend) = opbn(j)
349                        goto 20
350                     end if
351                  end do
352                  if (class(ib).eq.8 .or. class(ib).eq.17 .or.
353     &                class(ib).eq.26 .or. class(ib).eq.43 .or.
354     &                class(ib).eq.49 .or. class(ib).eq.73 .or.
355     &                class(ib).eq.82) then
356                     nopbend = nopbend + 1
357                     iopb(nopbend) = i
358                     opbk(nopbend) = 0.0d0
359                     goto 20
360                  end if
361                  goto 10
362   20             continue
363               end if
364            end if
365         end do
366      end if
367c
368c     turn off the out-of-plane bending term if it is not used
369c
370      if (nopbend .eq. 0)  use_opbend = .false.
371      return
372      end
373