wham and cluster_wham Adam's new constr_dist multichain
[unres.git] / source / unres / src_MD-M / 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 #ifdef DEBUG
275       write (iout,*) "Before sum_gradient"
276       do i=1,nres-1
277         write (iout,*) i," gradc  ",(gradc(j,i,icg),j=1,3)
278         write (iout,*) i," gradx  ",(gradx(j,i,icg),j=1,3)
279       enddo
280 #endif
281       call sum_gradient
282 #ifdef TIMING
283 #endif
284 #ifdef DEBUG
285       write (iout,*) "After sum_gradient"
286       do i=1,nres-1
287         write (iout,*) i," gradc  ",(gradc(j,i,icg),j=1,3)
288         write (iout,*) i," gradx  ",(gradx(j,i,icg),j=1,3)
289       enddo
290 #endif
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 #ifdef DEBUG
316       write (iout,*) "gcart, gxcart, gloc before int_to_cart"
317 #endif
318       do i=0,nct
319         do j=1,3
320           gcart(j,i)=gradc(j,i,icg)
321           gxcart(j,i)=gradx(j,i,icg)
322         enddo
323 #ifdef DEBUG
324         write (iout,'(i5,2(3f10.5,5x),f10.5)') i,(gcart(j,i),j=1,3),
325      &    (gxcart(j,i),j=1,3),gloc(i,icg)
326 #endif
327       enddo
328 #ifdef TIMING
329       time01=MPI_Wtime()
330 #endif
331       call int_to_cart
332 #ifdef TIMING
333       time_inttocart=time_inttocart+MPI_Wtime()-time01
334 #endif
335 #ifdef DEBUG
336       write (iout,*) "gcart and gxcart after int_to_cart"
337       do i=0,nres-1
338         write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
339      &      (gxcart(j,i),j=1,3)
340       enddo
341 #endif
342 #ifdef TIMING
343       time_cartgrad=time_cartgrad+MPI_Wtime()-time00
344 #endif
345       return
346       end
347 C-------------------------------------------------------------------------
348       subroutine zerograd
349       implicit real*8 (a-h,o-z)
350       include 'DIMENSIONS'
351       include 'COMMON.DERIV'
352       include 'COMMON.CHAIN'
353       include 'COMMON.VAR'
354       include 'COMMON.MD'
355       include 'COMMON.SCCOR'
356 C
357 C Initialize Cartesian-coordinate gradient
358 C
359       do i=-1,nres
360         do j=1,3
361           gvdwx(j,i)=0.0D0
362           gradx_scp(j,i)=0.0D0
363           gvdwc(j,i)=0.0D0
364           gvdwc_scp(j,i)=0.0D0
365           gvdwc_scpp(j,i)=0.0d0
366           gelc (j,i)=0.0D0
367           gelc_long(j,i)=0.0D0
368           gradb(j,i)=0.0d0
369           gradbx(j,i)=0.0d0
370           gvdwpp(j,i)=0.0d0
371           gel_loc(j,i)=0.0d0
372           gel_loc_long(j,i)=0.0d0
373           ghpbc(j,i)=0.0D0
374           ghpbx(j,i)=0.0D0
375           gsaxsc(j,i)=0.0D0
376           gsaxsx(j,i)=0.0D0
377           gcorr3_turn(j,i)=0.0d0
378           gcorr4_turn(j,i)=0.0d0
379           gradcorr(j,i)=0.0d0
380           gradcorr_long(j,i)=0.0d0
381           gradcorr5_long(j,i)=0.0d0
382           gradcorr6_long(j,i)=0.0d0
383           gcorr6_turn_long(j,i)=0.0d0
384           gradcorr5(j,i)=0.0d0
385           gradcorr6(j,i)=0.0d0
386           gcorr6_turn(j,i)=0.0d0
387           gsccorc(j,i)=0.0d0
388           gsccorx(j,i)=0.0d0
389           gradc(j,i,icg)=0.0d0
390           gradx(j,i,icg)=0.0d0
391           gscloc(j,i)=0.0d0
392           gsclocx(j,i)=0.0d0
393           gliptranc(j,i)=0.0d0
394           gliptranx(j,i)=0.0d0
395           gradafm(j,i)=0.0d0
396           do intertyp=1,3
397            gloc_sc(intertyp,i,icg)=0.0d0
398           enddo
399         enddo
400       enddo
401 c
402 c Initialize the gradients of local restraints
403 c
404       do i=1,nres
405         dutheta(i)=0.0d0
406         dugamma(i)=0.0d0
407         do j=1,3
408           duscdiff(j,i)=0.0d0
409           duscdiffx(j,i)=0.0d0
410         enddo
411       enddo
412 C
413 C Initialize the gradient of local energy terms.
414 C
415       do i=1,4*nres
416         gloc(i,icg)=0.0D0
417       enddo
418       do i=1,nres
419         gel_loc_loc(i)=0.0d0
420         gcorr_loc(i)=0.0d0
421         g_corr5_loc(i)=0.0d0
422         g_corr6_loc(i)=0.0d0
423         gel_loc_turn3(i)=0.0d0
424         gel_loc_turn4(i)=0.0d0
425         gel_loc_turn6(i)=0.0d0
426         gsccor_loc(i)=0.0d0
427       enddo
428 c initialize gcart and gxcart
429       do i=0,nres
430         do j=1,3
431           gcart(j,i)=0.0d0
432           gxcart(j,i)=0.0d0
433         enddo
434       enddo
435       return
436       end
437 c-------------------------------------------------------------------------
438       double precision function fdum()
439       fdum=0.0D0
440       return
441       end