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
71 gthetai=gthetai+dxdv(k,ind1)*gradx(k,j,icg)
72 gphii=gphii+dxdv(k+3,ind1)*gradx(k,j,icg)
75 if (i.gt.1) g(i-1)=gphii
76 if (n.gt.nphi) g(nphi+i)=gthetai
78 if (n.le.nphi+ntheta) goto 10
80 if (itype(i).ne.10) then
84 galphai=galphai+dxds(k,i)*gradx(k,i,icg)
87 gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
90 g(ialph(i,1)+nside)=gomegai
94 C Add the components corresponding to local energy terms.
97 c Add the usampl contributions
100 gloc(i,icg)=gloc(i,icg)+dugamma(i)
103 gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
107 cd write (iout,*) 'i=',i,'g=',g(i),' gloc=',gloc(i,icg)
108 g(i)=g(i)+gloc(i,icg)
110 C Uncomment following three lines for diagnostics.
112 cd call briefout(0,0.0d0)
113 cd write (iout,'(i3,1pe15.5)') (k,g(k),k=1,n)
116 C-------------------------------------------------------------------------
117 subroutine grad_restr(n,x,nf,g,uiparm,urparm,ufparm)
120 include 'COMMON.CHAIN'
121 include 'COMMON.DERIV'
123 include 'COMMON.INTERACT'
124 include 'COMMON.FFIELD'
125 include 'COMMON.IOUNITS'
127 double precision ufparm
130 double precision urparm(1)
131 double precision x(maxvar),g(maxvar)
132 integer i,j,k,ig,ind,ij,igall
133 double precision f,gthetai,gphii,galphai,gomegai
136 if (nf-nfl+1) 20,30,40
137 20 call func_restr(n,x,nf,f,uiparm,urparm,ufparm)
138 c write (iout,*) 'grad 20'
143 c Intercept NaNs in the coordinates
144 c write(iout,*) (var(i),i=1,nvar)
149 if (x_sum.ne.x_sum) then
150 write(iout,*)" *** grad_restr : Found NaN in coordinates"
152 print *," *** grad_restr : Found NaN in coordinates"
156 call var_to_geom_restr(n,x)
159 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
163 C Convert the Cartesian gradient into internal-coordinate gradient.
169 IF (mask_phi(i+2).eq.1) THEN
174 gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
175 gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
188 IF (mask_theta(i+2).eq.1) THEN
194 gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
195 gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
205 if (itype(i).ne.10) then
206 IF (mask_side(i).eq.1) THEN
210 galphai=galphai+dxds(k,i)*gradx(k,i,icg)
219 if (itype(i).ne.10) then
220 IF (mask_side(i).eq.1) THEN
224 gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
232 C Add the components corresponding to local energy terms.
239 if (mask_phi(i).eq.1) then
241 g(ig)=g(ig)+gloc(igall,icg)
247 if (mask_theta(i).eq.1) then
249 g(ig)=g(ig)+gloc(igall,icg)
255 if (itype(i).ne.10) then
257 if (mask_side(i).eq.1) then
259 g(ig)=g(ig)+gloc(igall,icg)
266 cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
270 C-------------------------------------------------------------------------
277 include 'COMMON.CONTROL'
278 include 'COMMON.CHAIN'
279 include 'COMMON.DERIV'
281 include 'COMMON.INTERACT'
282 include 'COMMON.FFIELD'
284 include 'COMMON.QRESTR'
285 include 'COMMON.IOUNITS'
286 include 'COMMON.TIME1'
289 c This subrouting calculates total Cartesian coordinate gradient.
290 c The subroutine chainbuild_cart and energy MUST be called beforehand.
297 write (iout,*) "Before sum_gradient"
299 write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3)
300 write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3)
302 write (iout,*) "gsaxsc, gsaxcx"
304 write (iout,*) i," gsaxsc ",(gsaxsc(j,i),j=1,3)
305 write (iout,*) i," gsaxsx ",(gsaxsx(j,i),j=1,3)
312 write (iout,*) "After sum_gradient"
314 write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3)
315 write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3)
318 c If performing constraint dynamics, add the gradients of the constraint energy
319 if(usampl.and.totT.gt.eq_time) then
321 write (iout,*) "dudconst, duscdiff, dugamma,dutheta"
322 write (iout,*) "wumb",wumb
324 write (iout,'(i5,3f10.5,5x,3f10.5,5x,2f10.5)')
325 & i,(dudconst(j,i),j=1,3),(duscdiff(j,i),j=1,3),
326 & dugamma(i),dutheta(i)
331 gradc(j,i,icg)=gradc(j,i,icg)+
332 & wumb*(dudconst(j,i)+duscdiff(j,i))
333 gradx(j,i,icg)=gradx(j,i,icg)+
334 & wumb*(dudxconst(j,i)+duscdiffx(j,i))
338 gloc(i,icg)=gloc(i,icg)+wumb*dugamma(i)
341 gloc(nphi+i,icg)=gloc(nphi+i,icg)+wumb*dutheta(i)
349 time_intcartderiv=time_intcartderiv+MPI_Wtime()-time01
351 cd call checkintcartgrad
352 cd write(iout,*) 'calling int_to_cart'
354 write (iout,*) "gcart, gxcart, gloc before int_to_cart"
358 gcart(j,i)=gradc(j,i,icg)
359 gxcart(j,i)=gradx(j,i,icg)
362 if((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
363 write (iout,'(i5,2(3f10.5,5x),4f10.5)') i,(gcart(j,i),j=1,3),
364 & (gxcart(j,i),j=1,3),gloc(i,icg),gloc(i+nphi,icg),
365 & gloc(ialph(i,1),icg),gloc(ialph(i,1)+nside,icg)
367 write (iout,'(i5,2(3f10.5,5x),4f10.5)') i,(gcart(j,i),j=1,3),
368 & (gxcart(j,i),j=1,3),gloc(i,icg),gloc(i+nphi,icg)
378 time_inttocart=time_inttocart+MPI_Wtime()-time01
381 write (iout,*) "gcart and gxcart after int_to_cart"
383 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
384 & (gxcart(j,i),j=1,3)
388 time_cartgrad=time_cartgrad+MPI_Wtime()-time00
392 c---------------------------------------------------------------------------
394 subroutine grad_transform
400 include 'COMMON.CONTROL'
401 include 'COMMON.CHAIN'
402 include 'COMMON.DERIV'
404 include 'COMMON.INTERACT'
405 include 'COMMON.FFIELD'
407 include 'COMMON.QRESTR'
408 include 'COMMON.IOUNITS'
409 include 'COMMON.TIME1'
412 write (iout,*)"Converting virtual-bond gradient to CA/SC gradient"
416 gcart(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
417 ! gcart_new(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
419 ! write (iout,'(i5,3f10.5,5x,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), &
420 ! (gcart_new(j,i),j=1,3),(gxcart(j,i),j=1,3)
422 ! Correction: dummy residues
425 gcart(j,nnt)=gcart(j,nnt)+gcart(j,1)
428 if (nct.lt.nres) then
430 ! gcart_new(j,nct)=gcart_new(j,nct)+gcart_new(j,nres)
431 gcart(j,nct)=gcart(j,nct)+gcart(j,nres)
435 write (iout,*) "CA/SC gradient"
437 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
438 & (gxcart(j,i),j=1,3)
444 C-------------------------------------------------------------------------
448 include 'COMMON.DERIV'
449 include 'COMMON.CHAIN'
452 include 'COMMON.SCCOR'
453 include 'COMMON.SHIELD'
454 integer i,j,kk,intertyp,maxshieldlist
457 C Initialize Cartesian-coordinate gradient
465 gvdwc_scpp(j,i)=0.0d0
467 C below is zero grad for shielding in order: ees (p-p)
468 C ecorr4, eturn3, eturn4, eel_loc, c denotes calfa,x is side-chain
471 gshieldc_loc(j,i)=0.0d0
472 gshieldx_ec(j,i)=0.0d0
473 gshieldc_ec(j,i)=0.0d0
474 gshieldc_loc_ec(j,i)=0.0d0
475 gshieldx_t3(j,i)=0.0d0
476 gshieldc_t3(j,i)=0.0d0
477 gshieldc_loc_t3(j,i)=0.0d0
478 gshieldx_t4(j,i)=0.0d0
479 gshieldc_t4(j,i)=0.0d0
480 gshieldc_loc_t4(j,i)=0.0d0
481 gshieldx_ll(j,i)=0.0d0
482 gshieldc_ll(j,i)=0.0d0
483 gshieldc_loc_ll(j,i)=0.0d0
484 C end of zero grad for shielding
490 gel_loc_long(j,i)=0.0d0
495 gcorr3_turn(j,i)=0.0d0
496 gcorr4_turn(j,i)=0.0d0
498 gradcorr_long(j,i)=0.0d0
499 gradcorr5_long(j,i)=0.0d0
500 gradcorr6_long(j,i)=0.0d0
501 gcorr6_turn_long(j,i)=0.0d0
504 gcorr6_turn(j,i)=0.0d0
514 grad_shield(j,i)=0.0d0
516 gg_tube_sc(j,i)=0.0d0
517 C grad_shield_side is Cbeta sidechain gradient
518 do kk=1,maxshieldlist
519 grad_shield_side(j,kk,i)=0.0d0
520 grad_shield_loc(j,kk,i)=0.0d0
522 C grad_shield_side_ca is Calfa sidechain gradient
525 C grad_shield_side_ca(j,kk,i)=0.0d0
528 gloc_sc(intertyp,i,icg)=0.0d0
543 C Initialize the gradient of local energy terms.
553 gel_loc_turn3(i)=0.0d0
554 gel_loc_turn4(i)=0.0d0
555 gel_loc_turn6(i)=0.0d0
558 c initialize gcart and gxcart
567 c-------------------------------------------------------------------------
568 double precision function fdum()