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