1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  2020  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     ################################################################
9c     ##                                                            ##
10c     ##  subroutine predict  --  induced dipole prediction values  ##
11c     ##                                                            ##
12c     ################################################################
13c
14c
15c     "predict" checks for use of methods for predicting induced
16c     dipoles, extrapolation coefficients and IELSCF parameters
17c
18c
19      subroutine predict
20      use atoms
21      use ielscf
22      use keys
23      use uprior
24      implicit none
25      integer i,j,k
26      integer next
27      character*20 keyword
28      character*240 record
29      character*240 string
30c
31c
32c     set defaults for use of induced dipole prediction
33c
34      use_pred = .false.
35      use_ielscf = .false.
36      polpred = '    '
37      maxualt = 0
38      nualt = 0
39c
40c     get keywords containing induced dipole prediction options
41c
42      do j = 1, nkey
43         next = 1
44         record = keyline(j)
45         call gettext (record,keyword,next)
46         call upcase (keyword)
47         string = record(next:240)
48         if (keyword(1:14) .eq. 'POLAR-PREDICT ') then
49            call getword (record,polpred,next)
50            call upcase (polpred)
51            use_pred = .true.
52            if (polpred .eq. '    ') then
53               polpred = 'ASPC'
54            else if (polpred .eq. 'IEL ') then
55               use_pred = .false.
56               use_ielscf = .true.
57            end if
58         else if (keyword(1:8) .eq. 'IEL-SCF ') then
59            use_ielscf = .true.
60         end if
61      end do
62c
63c     set always stable predictor-corrector (ASPC) coefficients
64c
65      if (polpred .eq. 'ASPC') then
66         maxualt = 17
67         aspc(1) = 62.0d0 / 17.0d0
68         aspc(2) = -310.0d0 / 51.0d0
69         aspc(3) = 2170.0d0 / 323.0d0
70         aspc(4) = -2329.0d0 / 400.0d0
71         aspc(5) = 1701.0d0 / 409.0d0
72         aspc(6) = -806.0d0 / 323.0d0
73         aspc(7) = 1024.0d0 / 809.0d0
74         aspc(8) = -479.0d0 / 883.0d0
75         aspc(9) = 257.0d0 / 1316.0d0
76         aspc(10) = -434.0d0 / 7429.0d0
77         aspc(11) = 191.0d0 / 13375.0d0
78         aspc(12) = -62.0d0 / 22287.0d0
79         aspc(13) = 3.0d0 / 7217.0d0
80         aspc(14) = -3.0d0 / 67015.0d0
81         aspc(15) = 2.0d0 / 646323.0d0
82         aspc(16) = -1.0d0 / 9694845.0d0
83         aspc(17) = 0.0d0
84      end if
85c
86c     set the 6th-order Gear predictor binomial coefficients
87c
88      if (polpred .eq. 'GEAR') then
89         maxualt = 7
90         gear(1) = 6.0d0
91         gear(2) = -15.0d0
92         gear(3) = 20.0d0
93         gear(4) = -15.0d0
94         gear(5) = 6.0d0
95         gear(6) = -1.0d0
96         gear(7) = 0.0d0
97      end if
98c
99c     set maximum storage size for least squares prediction
100c
101      if (polpred .eq. 'LSQR') then
102         maxualt = 6
103      end if
104c
105c     perform dynamic allocation of some global arrays
106c
107      if (use_pred) then
108         if (allocated(udalt))  deallocate (udalt)
109         if (allocated(upalt))  deallocate (upalt)
110         if (allocated(usalt))  deallocate (usalt)
111         if (allocated(upsalt))  deallocate (upsalt)
112         if (use_pred) then
113            allocate (udalt(maxualt,3,n))
114            allocate (upalt(maxualt,3,n))
115            allocate (usalt(maxualt,3,n))
116            allocate (upsalt(maxualt,3,n))
117         end if
118      end if
119c
120c     initialize prior values of induced dipole moments
121c
122      if (use_pred) then
123        do i = 1, n
124            do j = 1, 3
125               do k = 1, maxualt
126                  udalt(k,j,i) = 0.0d0
127                  upalt(k,j,i) = 0.0d0
128                  usalt(k,j,i) = 0.0d0
129                  upsalt(k,j,i) = 0.0d0
130               end do
131            end do
132         end do
133      end if
134c
135c     initialize inertial extended Lagrangian method
136c
137      if (use_ielscf)  call auxinit
138      return
139      end
140c
141c
142c     ##################################################################
143c     ##                                                              ##
144c     ##  subroutine auxinit  --  setup auxiliary dipoles for IELSCF  ##
145c     ##                                                              ##
146c     ##################################################################
147c
148c
149c     "auxinit" initializes auxiliary variables and settings for
150c     inertial extended Lagrangian induced dipole prediction
151c
152c     literature reference:
153c
154c     A. Albaugh, O. Demerdash, and T. Head-Gordon, "An Efficient and
155c     Stable Hybrid Extended Lagrangian/Self-Consistent Field Scheme
156c     for Solving Classical Mutual Induction", Journal of Chemical
157c     Physics, 143, 174104 (2015)
158c
159c
160      subroutine auxinit
161      use atomid
162      use atoms
163      use ielscf
164      use keys
165      use polar
166      implicit none
167      integer i,j,next
168      real*8 speed
169      real*8 weight
170      real*8 maxwell
171      real*8 vec(3)
172      character*20 keyword
173      character*240 record
174      character*240 string
175c
176c
177c     set defaults for auxiliary thermostat control variables
178c
179      nfree_aux = 3 * npolar
180      kelvin_aux = 100000.0d0
181      tautemp_aux = 0.1d0
182c
183c     check for keywords containing auxiliary thermostat values
184c
185      do i = 1, nkey
186         next = 1
187         record = keyline(i)
188         call gettext (record,keyword,next)
189         call upcase (keyword)
190         string = record(next:240)
191         if (keyword(1:13) .eq. 'AUX-TAUTEMP ') then
192            read (string,*,err=10,end=10)  tautemp_aux
193         else if (keyword(1:9) .eq. 'AUX-TEMP ') then
194            read (string,*,err=10,end=10)  kelvin_aux
195         end if
196   10    continue
197      end do
198c
199c     perform dynamic allocation of some global arrays
200c
201      allocate (uaux(3,n))
202      allocate (vaux(3,n))
203      allocate (aaux(3,n))
204      allocate (upaux(3,n))
205      allocate (vpaux(3,n))
206      allocate (apaux(3,n))
207c
208c     set auxiliary dipole values equal to induced dipoles
209c
210      use_ielscf = .false.
211      call induce
212      use_ielscf = .true.
213      do i = 1, n
214         do j = 1, 3
215            uaux(j,i) = uind(j,i)
216            upaux(j,i) = uinp(j,i)
217         end do
218      end do
219c
220c     set velocities and accelerations for auxiliary dipoles
221c
222      do i = 1, n
223         weight = 1.0d0
224         speed = maxwell (weight,kelvin_aux)
225         call ranvec (vec)
226         do j = 1, 3
227            vaux(j,i) = speed * vec(j)
228            aaux(j,i) = 0.0d0
229            vpaux(j,i) = vaux(j,i)
230            apaux(j,i) = 0.0d0
231         end do
232      end do
233      return
234      end
235