7fec1e8a0a1f3d2b75c3c91e4cc2924d49c0e5c0
[unres.git] / source / unres / src_MD / gradient_p.F
1       subroutine gradient(n,x,nf,g,uiparm,urparm,ufparm)
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4       include 'COMMON.CHAIN'
5       include 'COMMON.DERIV'
6       include 'COMMON.VAR'
7       include 'COMMON.INTERACT'
8       include 'COMMON.FFIELD'
9       include 'COMMON.MD'
10       include 'COMMON.IOUNITS'
11       include 'COMMON.SCCOR'
12       external ufparm
13       integer uiparm(1)
14       double precision urparm(1)
15       dimension x(maxvar),g(maxvar)
16 c
17 c This subroutine calculates total internal coordinate gradient.
18 c Depending on the number of function evaluations, either whole energy 
19 c is evaluated beforehand, Cartesian coordinates and their derivatives in 
20 c internal coordinates are reevaluated or only the cartesian-in-internal
21 c coordinate derivatives are evaluated. The subroutine was designed to work
22 c with SUMSL.
23
24 c
25       icg=mod(nf,2)+1
26
27 cd      print *,'grad',nf,icg
28       if (nf-nfl+1) 20,30,40
29    20 call func(n,x,nf,f,uiparm,urparm,ufparm)
30 c     write (iout,*) 'grad 20'
31       if (nf.eq.0) return
32       goto 40
33    30 call var_to_geom(n,x)
34       call chainbuild 
35 c     write (iout,*) 'grad 30'
36 C
37 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
38 C
39    40 call cartder
40 c     write (iout,*) 'grad 40'
41 c     print *,'GRADIENT: nnt=',nnt,' nct=',nct,' expon=',expon
42 C
43 C Convert the Cartesian gradient into internal-coordinate gradient.
44 C
45       ind=0
46       ind1=0
47       do i=1,nres-2
48         gthetai=0.0D0
49         gphii=0.0D0
50         do j=i+1,nres-1
51           ind=ind+1
52 c         ind=indmat(i,j)
53 c         print *,'GRAD: i=',i,' jc=',j,' ind=',ind
54           do k=1,3
55             gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
56           enddo
57           do k=1,3
58             gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
59           enddo
60         enddo
61         do j=i+1,nres-1
62           ind1=ind1+1
63 c         ind1=indmat(i,j)
64 c         print *,'GRAD: i=',i,' jx=',j,' ind1=',ind1
65           do k=1,3
66             gthetai=gthetai+dxdv(k,ind1)*gradx(k,j,icg)
67             gphii=gphii+dxdv(k+3,ind1)*gradx(k,j,icg)
68           enddo
69         enddo
70         if (i.gt.1) g(i-1)=gphii
71         if (n.gt.nphi) g(nphi+i)=gthetai
72       enddo
73       if (n.le.nphi+ntheta) goto 10
74       do i=2,nres-1
75         if (itype(i).ne.10) then
76           galphai=0.0D0
77           gomegai=0.0D0
78           do k=1,3
79             galphai=galphai+dxds(k,i)*gradx(k,i,icg)
80           enddo
81           do k=1,3
82             gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
83           enddo
84           g(ialph(i,1))=galphai
85           g(ialph(i,1)+nside)=gomegai
86         endif
87       enddo
88 C
89 C Add the components corresponding to local energy terms.
90 C
91    10 continue
92       do i=1,nvar
93 cd      write (iout,*) 'i=',i,'g=',g(i),' gloc=',gloc(i,icg)
94         g(i)=g(i)+gloc(i,icg)
95       enddo
96 C Uncomment following three lines for diagnostics.
97 cd    call intout
98 cd    call briefout(0,0.0d0)
99 cd    write (iout,'(i3,1pe15.5)') (k,g(k),k=1,n)
100       return
101       end
102 C-------------------------------------------------------------------------
103       subroutine grad_restr(n,x,nf,g,uiparm,urparm,ufparm)
104       implicit real*8 (a-h,o-z)
105       include 'DIMENSIONS'
106       include 'COMMON.CHAIN'
107       include 'COMMON.DERIV'
108       include 'COMMON.VAR'
109       include 'COMMON.INTERACT'
110       include 'COMMON.FFIELD'
111       include 'COMMON.IOUNITS'
112       external ufparm
113       integer uiparm(1)
114       double precision urparm(1)
115       dimension x(maxvar),g(maxvar)
116
117       icg=mod(nf,2)+1
118       if (nf-nfl+1) 20,30,40
119    20 call func_restr(n,x,nf,f,uiparm,urparm,ufparm)
120 c     write (iout,*) 'grad 20'
121       if (nf.eq.0) return
122       goto 40
123    30 continue
124 #ifdef OSF
125 c     Intercept NaNs in the coordinates
126 c      write(iout,*) (var(i),i=1,nvar)
127       x_sum=0.D0
128       do i=1,n
129         x_sum=x_sum+x(i)
130       enddo
131       if (x_sum.ne.x_sum) then
132         write(iout,*)" *** grad_restr : Found NaN in coordinates"
133         call flush(iout)
134         print *," *** grad_restr : Found NaN in coordinates"
135         return
136       endif
137 #endif
138       call var_to_geom_restr(n,x)
139       call chainbuild 
140 C
141 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
142 C
143    40 call cartder
144 C
145 C Convert the Cartesian gradient into internal-coordinate gradient.
146 C
147
148       ig=0
149       ind=nres-2                                                                    
150       do i=2,nres-2                
151        IF (mask_phi(i+2).eq.1) THEN                                             
152         gphii=0.0D0                                                             
153         do j=i+1,nres-1                                                         
154           ind=ind+1                                 
155           do k=1,3                                                              
156             gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)                            
157             gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)                           
158           enddo                                                                 
159         enddo                                                                   
160         ig=ig+1
161         g(ig)=gphii
162        ELSE
163         ind=ind+nres-1-i
164        ENDIF
165       enddo                                        
166
167
168       ind=0
169       do i=1,nres-2
170        IF (mask_theta(i+2).eq.1) THEN
171         ig=ig+1
172         gthetai=0.0D0
173         do j=i+1,nres-1
174           ind=ind+1
175           do k=1,3
176             gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
177             gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
178           enddo
179         enddo
180         g(ig)=gthetai
181        ELSE
182         ind=ind+nres-1-i
183        ENDIF
184       enddo
185
186       do i=2,nres-1
187         if (itype(i).ne.10) then
188          IF (mask_side(i).eq.1) THEN
189           ig=ig+1
190           galphai=0.0D0
191           do k=1,3
192             galphai=galphai+dxds(k,i)*gradx(k,i,icg)
193           enddo
194           g(ig)=galphai
195          ENDIF
196         endif
197       enddo
198
199       
200       do i=2,nres-1
201         if (itype(i).ne.10) then
202          IF (mask_side(i).eq.1) THEN
203           ig=ig+1
204           gomegai=0.0D0
205           do k=1,3
206             gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
207           enddo
208           g(ig)=gomegai
209          ENDIF
210         endif
211       enddo
212
213 C
214 C Add the components corresponding to local energy terms.
215 C
216
217       ig=0
218       igall=0
219       do i=4,nres
220         igall=igall+1
221         if (mask_phi(i).eq.1) then
222           ig=ig+1
223           g(ig)=g(ig)+gloc(igall,icg)
224         endif
225       enddo
226
227       do i=3,nres
228         igall=igall+1
229         if (mask_theta(i).eq.1) then
230           ig=ig+1
231           g(ig)=g(ig)+gloc(igall,icg)
232         endif
233       enddo
234      
235       do ij=1,2
236       do i=2,nres-1
237         if (itype(i).ne.10) then
238           igall=igall+1
239           if (mask_side(i).eq.1) then
240             ig=ig+1
241             g(ig)=g(ig)+gloc(igall,icg)
242           endif
243         endif
244       enddo
245       enddo
246
247 cd      do i=1,ig
248 cd        write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
249 cd      enddo
250       return
251       end
252 C-------------------------------------------------------------------------
253       subroutine cartgrad
254       implicit real*8 (a-h,o-z)
255       include 'DIMENSIONS'
256 #ifdef MPI
257       include 'mpif.h'
258 #endif
259       include 'COMMON.CHAIN'
260       include 'COMMON.DERIV'
261       include 'COMMON.VAR'
262       include 'COMMON.INTERACT'
263       include 'COMMON.FFIELD'
264       include 'COMMON.MD'
265       include 'COMMON.IOUNITS'
266       include 'COMMON.TIME1'
267       include 'COMMON.SCCOR'
268 c
269 c This subrouting calculates total Cartesian coordinate gradient. 
270 c The subroutine chainbuild_cart and energy MUST be called beforehand.
271 c
272 c        do i=1,nres
273 c        write (iout,*) "przed sum_grad", gloc_sc(1,i,icg),gloc(i,icg)
274 c        enddo
275
276 #ifdef TIMING
277       time00=MPI_Wtime()
278 #endif
279       icg=1
280       call sum_gradient
281 #ifdef TIMING
282 #endif
283 c        do i=1,nres
284 c        write (iout,*) "checkgrad", gloc_sc(1,i,icg),gloc(i,icg)
285 c        enddo     
286 cd      write (iout,*) "After sum_gradient"
287 cd      do i=1,nres-1
288 cd        write (iout,*) i," gradc  ",(gradc(j,i,icg),j=1,3)
289 cd        write (iout,*) i," gradx  ",(gradx(j,i,icg),j=1,3)
290 cd      enddo
291 c If performing constraint dynamics, add the gradients of the constraint energy
292       if(usampl.and.totT.gt.eq_time) then
293          do i=1,nct
294            do j=1,3
295              gradc(j,i,icg)=gradc(j,i,icg)+dudconst(j,i)+duscdiff(j,i)
296              gradx(j,i,icg)=gradx(j,i,icg)+dudxconst(j,i)+duscdiffx(j,i)
297            enddo
298          enddo
299          do i=1,nres-3
300            gloc(i,icg)=gloc(i,icg)+dugamma(i)
301          enddo
302          do i=1,nres-2
303            gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
304          enddo
305       endif 
306 #ifdef TIMING
307       time01=MPI_Wtime()
308 #endif
309       call intcartderiv
310 #ifdef TIMING
311       time_intcartderiv=time_intcartderiv+MPI_Wtime()-time01
312 #endif
313 cd      call checkintcartgrad
314 cd      write(iout,*) 'calling int_to_cart'
315 cd      write (iout,*) "gcart, gxcart, gloc before int_to_cart"
316       do i=1,nct
317         do j=1,3
318           gcart(j,i)=gradc(j,i,icg)
319           gxcart(j,i)=gradx(j,i,icg)
320         enddo
321 cd        write (iout,'(i5,2(3f10.5,5x),f10.5)') i,(gcart(j,i),j=1,3),
322 cd     &    (gxcart(j,i),j=1,3),gloc(i,icg)
323       enddo
324 #ifdef TIMING
325       time01=MPI_Wtime()
326 #endif
327       call int_to_cart
328 #ifdef TIMING
329       time_inttocart=time_inttocart+MPI_Wtime()-time01
330 #endif
331 cd      write (iout,*) "gcart and gxcart after int_to_cart"
332 cd      do i=0,nres-1
333 cd        write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
334 cd     &      (gxcart(j,i),j=1,3)
335 cd      enddo
336 #ifdef TIMING
337       time_cartgrad=time_cartgrad+MPI_Wtime()-time00
338 #endif
339       return
340       end
341 C-------------------------------------------------------------------------
342       subroutine zerograd
343       implicit real*8 (a-h,o-z)
344       include 'DIMENSIONS'
345       include 'COMMON.DERIV'
346       include 'COMMON.CHAIN'
347       include 'COMMON.VAR'
348       include 'COMMON.MD'
349       include 'COMMON.SCCOR'
350 C
351 C Initialize Cartesian-coordinate gradient
352 C
353       do i=1,nres
354         do j=1,3
355           gvdwx(j,i)=0.0D0
356           gvdwxT(j,i)=0.0D0
357           gradx_scp(j,i)=0.0D0
358           gvdwc(j,i)=0.0D0
359           gvdwcT(j,i)=0.0D0
360           gvdwc_scp(j,i)=0.0D0
361           gvdwc_scpp(j,i)=0.0d0
362           gelc (j,i)=0.0D0
363           gelc_long(j,i)=0.0D0
364           gradb(j,i)=0.0d0
365           gradbx(j,i)=0.0d0
366           gvdwpp(j,i)=0.0d0
367           gel_loc(j,i)=0.0d0
368           gel_loc_long(j,i)=0.0d0
369           ghpbc(j,i)=0.0D0
370           ghpbx(j,i)=0.0D0
371           gcorr3_turn(j,i)=0.0d0
372           gcorr4_turn(j,i)=0.0d0
373           gradcorr(j,i)=0.0d0
374           gradcorr_long(j,i)=0.0d0
375           gradcorr5_long(j,i)=0.0d0
376           gradcorr6_long(j,i)=0.0d0
377           gcorr6_turn_long(j,i)=0.0d0
378           gradcorr5(j,i)=0.0d0
379           gradcorr6(j,i)=0.0d0
380           gcorr6_turn(j,i)=0.0d0
381           gsccorc(j,i)=0.0d0
382           gsccorx(j,i)=0.0d0
383           gradc(j,i,icg)=0.0d0
384           gradx(j,i,icg)=0.0d0
385           gscloc(j,i)=0.0d0
386           gsclocx(j,i)=0.0d0
387           do intertyp=1,3
388            gloc_sc(intertyp,i,icg)=0.0d0
389           enddo
390         enddo
391       enddo
392 C
393 C Initialize the gradient of local energy terms.
394 C
395       do i=1,4*nres
396         gloc(i,icg)=0.0D0
397       enddo
398       do i=1,nres
399         gel_loc_loc(i)=0.0d0
400         gcorr_loc(i)=0.0d0
401         g_corr5_loc(i)=0.0d0
402         g_corr6_loc(i)=0.0d0
403         gel_loc_turn3(i)=0.0d0
404         gel_loc_turn4(i)=0.0d0
405         gel_loc_turn6(i)=0.0d0
406         gsccor_loc(i)=0.0d0
407       enddo
408 c initialize gcart and gxcart
409       do i=0,nres
410         do j=1,3
411           gcart(j,i)=0.0d0
412           gxcart(j,i)=0.0d0
413         enddo
414       enddo
415       return
416       end
417 c-------------------------------------------------------------------------
418       double precision function fdum()
419       fdum=0.0D0
420       return
421       end