1 subroutine gradient(n,x,nf,g,uiparm,urparm,ufparm)
4 include 'COMMON.CONTROL'
8 include 'COMMON.INTERACT'
9 include 'COMMON.FFIELD'
11 include 'COMMON.QRESTR'
12 include 'COMMON.IOUNITS'
14 double precision ufparm
17 double precision urparm(1)
18 double precision x(n),g(n)
19 integer i,j,k,ind,ind1
20 double precision f,gthetai,gphii,galphai,gomegai
22 c This subroutine calculates total internal coordinate gradient.
23 c Depending on the number of function evaluations, either whole energy
24 c is evaluated beforehand, Cartesian coordinates and their derivatives in
25 c internal coordinates are reevaluated or only the cartesian-in-internal
26 c coordinate derivatives are evaluated. The subroutine was designed to work
32 cd print *,'grad',nf,icg
33 if (nf-nfl+1) 20,30,40
34 20 call func(n,x,nf,f,uiparm,urparm,ufparm)
35 c write (iout,*) 'grad 20'
38 30 call var_to_geom(n,x)
40 c write (iout,*) 'grad 30'
42 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
45 c write (iout,*) 'grad 40'
46 c print *,'GRADIENT: nnt=',nnt,' nct=',nct,' expon=',expon
48 C Convert the Cartesian gradient into internal-coordinate gradient.
58 c print *,'GRAD: i=',i,' jc=',j,' ind=',ind
60 gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
63 gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
69 c print *,'GRAD: i=',i,' jx=',j,' ind1=',ind1
70 write (iout,*) "i",i," j",j," ind1",ind1
71 write (iout,*) "dxdv",(dxdv(k,ind1),k=1,6)
72 write (iout,*) "gradx",(gradx(k,j,icg),k=1,3)
74 gthetai=gthetai+dxdv(k,ind1)*gradx(k,j,icg)
75 gphii=gphii+dxdv(k+3,ind1)*gradx(k,j,icg)
78 if (i.gt.1) g(i-1)=gphii
79 if (n.gt.nphi) g(nphi+i)=gthetai
81 if (n.le.nphi+ntheta) goto 10
83 if (itype(i).ne.10) then
87 galphai=galphai+dxds(k,i)*gradx(k,i,icg)
90 gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
93 g(ialph(i,1)+nside)=gomegai
97 C Add the components corresponding to local energy terms.
100 c Add the usampl contributions
103 gloc(i,icg)=gloc(i,icg)+dugamma(i)
106 gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
110 cd write (iout,*) 'i=',i,'g=',g(i),' gloc=',gloc(i,icg)
111 g(i)=g(i)+gloc(i,icg)
113 C Uncomment following three lines for diagnostics.
115 cd call briefout(0,0.0d0)
116 cd write (iout,'(i3,1pe15.5)') (k,g(k),k=1,n)
119 C-------------------------------------------------------------------------
120 subroutine grad_restr(n,x,nf,g,uiparm,urparm,ufparm)
123 include 'COMMON.CHAIN'
124 include 'COMMON.DERIV'
126 include 'COMMON.INTERACT'
127 include 'COMMON.FFIELD'
128 include 'COMMON.IOUNITS'
130 double precision ufparm
133 double precision urparm(1)
134 double precision x(maxvar),g(maxvar)
135 integer i,j,k,ig,ind,ij,igall
136 double precision f,gthetai,gphii,galphai,gomegai
139 if (nf-nfl+1) 20,30,40
140 20 call func_restr(n,x,nf,f,uiparm,urparm,ufparm)
141 c write (iout,*) 'grad 20'
146 c Intercept NaNs in the coordinates
147 c write(iout,*) (var(i),i=1,nvar)
152 if (x_sum.ne.x_sum) then
153 write(iout,*)" *** grad_restr : Found NaN in coordinates"
155 print *," *** grad_restr : Found NaN in coordinates"
159 call var_to_geom_restr(n,x)
162 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
166 C Convert the Cartesian gradient into internal-coordinate gradient.
172 IF (mask_phi(i+2).eq.1) THEN
177 gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
178 gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
191 IF (mask_theta(i+2).eq.1) THEN
197 gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
198 gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
208 if (itype(i).ne.10) then
209 IF (mask_side(i).eq.1) THEN
213 galphai=galphai+dxds(k,i)*gradx(k,i,icg)
222 if (itype(i).ne.10) then
223 IF (mask_side(i).eq.1) THEN
227 gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
235 C Add the components corresponding to local energy terms.
242 if (mask_phi(i).eq.1) then
244 g(ig)=g(ig)+gloc(igall,icg)
250 if (mask_theta(i).eq.1) then
252 g(ig)=g(ig)+gloc(igall,icg)
258 if (itype(i).ne.10) then
260 if (mask_side(i).eq.1) then
262 g(ig)=g(ig)+gloc(igall,icg)
269 cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
273 C-------------------------------------------------------------------------
280 include 'COMMON.CONTROL'
281 include 'COMMON.CHAIN'
282 include 'COMMON.DERIV'
284 include 'COMMON.INTERACT'
285 include 'COMMON.FFIELD'
287 include 'COMMON.QRESTR'
288 include 'COMMON.IOUNITS'
289 include 'COMMON.TIME1'
292 c This subrouting calculates total Cartesian coordinate gradient.
293 c The subroutine chainbuild_cart and energy MUST be called beforehand.
300 write (iout,*) "Before sum_gradient"
302 write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3)
303 write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3)
305 write (iout,*) "gsaxsc, gsaxcx"
307 write (iout,*) i," gsaxsc ",(gsaxsc(j,i),j=1,3)
308 write (iout,*) i," gsaxsx ",(gsaxsx(j,i),j=1,3)
315 write (iout,*) "After sum_gradient"
317 write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3)
318 write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3)
321 c If performing constraint dynamics, add the gradients of the constraint energy
322 if(usampl.and.totT.gt.eq_time) then
324 write (iout,*) "dudconst, duscdiff, dugamma,dutheta"
325 write (iout,*) "wumb",wumb
327 write (iout,'(i5,3f10.5,5x,3f10.5,5x,2f10.5)')
328 & i,(dudconst(j,i),j=1,3),(duscdiff(j,i),j=1,3),
329 & dugamma(i),dutheta(i)
334 gradc(j,i,icg)=gradc(j,i,icg)+
335 & wumb*(dudconst(j,i)+duscdiff(j,i))
336 gradx(j,i,icg)=gradx(j,i,icg)+
337 & wumb*(dudxconst(j,i)+duscdiffx(j,i))
341 gloc(i,icg)=gloc(i,icg)+wumb*dugamma(i)
344 gloc(nphi+i,icg)=gloc(nphi+i,icg)+wumb*dutheta(i)
352 time_intcartderiv=time_intcartderiv+MPI_Wtime()-time01
354 cd call checkintcartgrad
355 cd write(iout,*) 'calling int_to_cart'
357 write (iout,*) "gcart, gxcart, gloc before int_to_cart"
361 gcart(j,i)=gradc(j,i,icg)
362 gxcart(j,i)=gradx(j,i,icg)
365 if((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
366 write (iout,'(i5,2(3f10.5,5x),4f10.5)') i,(gcart(j,i),j=1,3),
367 & (gxcart(j,i),j=1,3),gloc(i,icg),gloc(i+nphi,icg),
368 & gloc(ialph(i,1),icg),gloc(ialph(i,1)+nside,icg)
370 write (iout,'(i5,2(3f10.5,5x),4f10.5)') i,(gcart(j,i),j=1,3),
371 & (gxcart(j,i),j=1,3),gloc(i,icg),gloc(i+nphi,icg)
381 time_inttocart=time_inttocart+MPI_Wtime()-time01
384 write (iout,*) "gcart and gxcart after int_to_cart"
386 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
387 & (gxcart(j,i),j=1,3)
391 time_cartgrad=time_cartgrad+MPI_Wtime()-time00
395 c---------------------------------------------------------------------------
397 subroutine grad_transform
403 include 'COMMON.CONTROL'
404 include 'COMMON.CHAIN'
405 include 'COMMON.DERIV'
407 include 'COMMON.INTERACT'
408 include 'COMMON.FFIELD'
410 include 'COMMON.QRESTR'
411 include 'COMMON.IOUNITS'
412 include 'COMMON.TIME1'
415 write (iout,*)"Converting virtual-bond gradient to CA/SC gradient"
419 gcart(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
420 ! gcart_new(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
422 ! write (iout,'(i5,3f10.5,5x,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), &
423 ! (gcart_new(j,i),j=1,3),(gxcart(j,i),j=1,3)
425 ! Correction: dummy residues
428 gcart(j,nnt)=gcart(j,nnt)+gcart(j,1)
431 if (nct.lt.nres) then
433 ! gcart_new(j,nct)=gcart_new(j,nct)+gcart_new(j,nres)
434 gcart(j,nct)=gcart(j,nct)+gcart(j,nres)
438 write (iout,*) "CA/SC gradient"
440 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
441 & (gxcart(j,i),j=1,3)
447 C-------------------------------------------------------------------------
451 include 'COMMON.DERIV'
452 include 'COMMON.CHAIN'
455 include 'COMMON.SCCOR'
456 include 'COMMON.SHIELD'
457 integer i,j,kk,intertyp,maxshieldlist
460 C Initialize Cartesian-coordinate gradient
468 gvdwc_scpp(j,i)=0.0d0
470 C below is zero grad for shielding in order: ees (p-p)
471 C ecorr4, eturn3, eturn4, eel_loc, c denotes calfa,x is side-chain
474 gshieldc_loc(j,i)=0.0d0
475 gshieldx_ec(j,i)=0.0d0
476 gshieldc_ec(j,i)=0.0d0
477 gshieldc_loc_ec(j,i)=0.0d0
478 gshieldx_t3(j,i)=0.0d0
479 gshieldc_t3(j,i)=0.0d0
480 gshieldc_loc_t3(j,i)=0.0d0
481 gshieldx_t4(j,i)=0.0d0
482 gshieldc_t4(j,i)=0.0d0
483 gshieldc_loc_t4(j,i)=0.0d0
484 gshieldx_ll(j,i)=0.0d0
485 gshieldc_ll(j,i)=0.0d0
486 gshieldc_loc_ll(j,i)=0.0d0
487 C end of zero grad for shielding
493 gel_loc_long(j,i)=0.0d0
498 gcorr3_turn(j,i)=0.0d0
499 gcorr4_turn(j,i)=0.0d0
501 gradcorr_long(j,i)=0.0d0
502 gradcorr5_long(j,i)=0.0d0
503 gradcorr6_long(j,i)=0.0d0
504 gcorr6_turn_long(j,i)=0.0d0
507 gcorr6_turn(j,i)=0.0d0
517 grad_shield(j,i)=0.0d0
519 gg_tube_sc(j,i)=0.0d0
520 C grad_shield_side is Cbeta sidechain gradient
521 do kk=1,maxshieldlist
522 grad_shield_side(j,kk,i)=0.0d0
523 grad_shield_loc(j,kk,i)=0.0d0
525 C grad_shield_side_ca is Calfa sidechain gradient
528 C grad_shield_side_ca(j,kk,i)=0.0d0
531 gloc_sc(intertyp,i,icg)=0.0d0
546 C Initialize the gradient of local energy terms.
556 gel_loc_turn3(i)=0.0d0
557 gel_loc_turn4(i)=0.0d0
558 gel_loc_turn6(i)=0.0d0
561 c initialize gcart and gxcart
570 c-------------------------------------------------------------------------
571 double precision function fdum()