1 C----------------------------------------------------------------------------
3 C Check the gradient of the energy in Cartesian coordinates.
6 include 'COMMON.CONTROL'
9 include 'COMMON.IOUNITS'
11 include 'COMMON.CONTACTS'
15 double precision ggg(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),
16 & g(maxvar),grad_s(6,maxres)
17 double precision energia(0:n_ene),energia1(0:n_ene)
18 double precision aincr2,etot,etot1,etot2
19 double precision dist,alpha,beta
20 double precision funcgrad,ff
24 double precision urparm(1)
31 print '("Calling CHECK_ECART",1pd12.3)',aincr
32 write (iout,'("Calling CHECK_ECART",1pd12.3)') aincr
35 call geom_to_var(nvar,x)
36 call etotal(energia(0))
38 call enerprint(energia(0))
42 call gradient(nvar,x,nf,g,uiparm,urparm,fdum)
46 write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
50 grad_s(j,i)=gradc(j,i,icg)
51 grad_s(j+3,i)=gradx(j,i,icg)
55 write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
66 c(j,k+nres)=c(j,k+nres)+aincr
68 call etotal(energia1(0))
70 ggg(j)=(etot1-etot)/aincr
74 c(j,k+nres)=c(j,k+nres)-aincr
78 c(j,i+nres)=c(j,i+nres)+aincr
79 dc(j,i+nres)=dc(j,i+nres)+aincr
80 call etotal(energia1(0))
82 ggg(j+3)=(etot1-etot)/aincr
86 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/)')
87 & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6)
91 c----------------------------------------------------------------------------
92 subroutine check_ecartint
93 C Check the gradient of the energy in Cartesian coordinates.
96 include 'COMMON.CONTROL'
97 include 'COMMON.CHAIN'
98 include 'COMMON.DERIV'
99 include 'COMMON.IOUNITS'
101 include 'COMMON.CONTACTS'
103 include 'COMMON.LOCAL'
104 include 'COMMON.SPLITELE'
107 double precision ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3),
108 & x(maxvar),g(maxvar)
109 double precision dcnorm_safe(3),dxnorm_safe(3)
110 double precision grad_s(6,0:maxres),grad_s1(6,0:maxres)
111 double precision phi_temp(maxres),theta_temp(maxres),
112 & alph_temp(maxres),omeg_temp(maxres)
113 double precision energia(0:n_ene),energia1(0:n_ene)
115 double precision urparm(1)
118 double precision etot,etot1,etot2,etot11,etot12,etot21,etot22
119 double precision dist,alpha,beta
126 call int_from_cart1(.false.)
129 c call checkintcartgrad
133 print '("Calling CHECK_ECARTINT",1pd12.3)',aincr
134 write (iout,'("Calling CHECK_ECARTINT",1pd12.3)') aincr
137 call geom_to_var(nvar,x)
138 if (.not.split_ene) then
139 call etotal(energia(0))
141 call enerprint(energia(0))
142 c write (iout,*) "enter cartgrad"
145 c write (iout,*) "exit cartgrad"
148 write (iout,'(//27(1h*)," Checking energy gradient ",27(1h*))')
149 write (iout,'(//4x,3a12,3x,3a12)')"gcart_x","gcart_y","gcart_z",
150 & "gxcart_x","gxcart_y","gxcart_z"
152 write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
153 & (gxcart(j,i),j=1,3)
156 grad_s(j,0)=gcart(j,0)
160 grad_s(j,i)=gcart(j,i)
161 grad_s(j+3,i)=gxcart(j,i)
165 !- split gradient check
167 call etotal_long(energia(0))
168 call enerprint(energia(0))
170 c write (iout,*) "enter cartgrad"
173 c write (iout,*) "exit cartgrad"
176 write (iout,*) "longrange grad"
178 write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
179 & (gxcart(j,i),j=1,3)
182 grad_s(j,0)=gcart(j,0)
186 grad_s(j,i)=gcart(j,i)
187 grad_s(j+3,i)=gxcart(j,i)
191 call etotal_short(energia(0))
192 call enerprint(energia(0))
194 c write (iout,*) "enter cartgrad"
197 c write (iout,*) "exit cartgrad"
200 write (iout,*) "shortrange grad"
202 write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
203 & (gxcart(j,i),j=1,3)
206 grad_s1(j,0)=gcart(j,0)
210 grad_s1(j,i)=gcart(j,i)
211 grad_s1(j+3,i)=gxcart(j,i)
215 write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
223 dcnorm_safe(k)=dc_norm(k,i)
224 dxnorm_safe(k)=dc_norm(k,i+nres)
231 c Broadcast the order to compute internal coordinates to the slaves.
233 c & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
235 c call int_from_cart1(.false.)
236 if (.not.split_ene) then
237 call etotal(energia1(0))
241 call etotal_long(energia1(0))
243 call etotal_short(energia1(0))
245 c write (iout,*) "etot11",etot11," etot12",etot12
247 !- end split gradient
248 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
252 c call int_from_cart1(.false.)
253 if (.not.split_ene) then
254 call etotal(energia1(0))
256 ggg(j)=(etot1-etot2)/(2*aincr)
259 call etotal_long(energia1(0))
261 ggg(j)=(etot11-etot21)/(2*aincr)
262 call etotal_short(energia1(0))
264 ggg1(j)=(etot12-etot22)/(2*aincr)
265 !- end split gradient
266 c write (iout,*) "etot21",etot21," etot22",etot22
268 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
273 dc(j,i+nres)=ddx(j)+aincr
275 c write (iout,*) "i",i," j",j," dxnorm+ and dxnorm"
276 c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
277 c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
278 c write (iout,*) "dxnormnorm",dsqrt(
279 c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
280 c write (iout,*) "dxnormnormsafe",dsqrt(
281 c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
283 if (.not.split_ene) then
284 call etotal(energia1(0))
288 call etotal_long(energia1(0))
290 call etotal_short(energia1(0))
293 !- end split gradient
294 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
295 dc(j,i+nres)=ddx(j)-aincr
297 c write (iout,*) "i",i," j",j," dxnorm- and dxnorm"
298 c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
299 c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
301 c write (iout,*) "dxnormnorm",dsqrt(
302 c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
303 c write (iout,*) "dxnormnormsafe",dsqrt(
304 c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
305 if (.not.split_ene) then
306 call etotal(energia1(0))
308 ggg(j+3)=(etot1-etot2)/(2*aincr)
311 call etotal_long(energia1(0))
313 ggg(j+3)=(etot11-etot21)/(2*aincr)
314 call etotal_short(energia1(0))
316 ggg1(j+3)=(etot12-etot22)/(2*aincr)
317 !- end split gradient
319 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
323 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
324 & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
326 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
327 & i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),
329 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
330 & i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),
331 & ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
336 c-------------------------------------------------------------------------
337 subroutine int_from_cart1(lprn)
344 include 'COMMON.IOUNITS'
346 include 'COMMON.CHAIN'
348 include 'COMMON.INTERACT'
349 include 'COMMON.LOCAL'
350 include 'COMMON.NAMES'
351 include 'COMMON.SETUP'
352 include 'COMMON.TIME1'
355 double precision dnorm1,dnorm2,be
356 double precision time00
357 double precision dist,alpha,beta
358 if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
362 #if defined(PARINT) && defined(MPI)
363 do i=iint_start,iint_end
370 C print *,i,dnorm1,dnorm2
372 c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1
373 & +(c(j,i+1)-c(j,i))/dnorm2)
377 if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
378 if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
379 tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
381 if (itype(i-1).ne.10) then
382 tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
383 omicron(1,i)=alpha(i-2,i-1,i-1+nres)
384 omicron(2,i)=alpha(i-1+nres,i-1,i)
386 if (itype(i).ne.10) then
387 tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
390 omeg(i)=beta(nres+i,i,maxres2,i+1)
392 alph(i)=alpha(nres+i,i,maxres2)
394 theta(i+1)=alpha(i-1,i,i+1)
397 vbld_inv(i)=1.0d0/vbld(i)
398 vbld(nres+i)=dist(nres+i,i)
399 C print *,vbld(i+nres)
401 if (itype(i).ne.10) then
402 vbld_inv(nres+i)=1.0d0/vbld(nres+i)
404 vbld_inv(nres+i)=0.0d0
407 #if defined(PARINT) && defined(MPI)
408 if (nfgtasks1.gt.1) then
409 cd write(iout,*) "iint_start",iint_start," iint_count",
410 cd & (iint_count(i),i=0,nfgtasks-1)," iint_displ",
411 cd & (iint_displ(i),i=0,nfgtasks-1)
412 cd write (iout,*) "Gather vbld backbone"
415 call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),
416 & MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),
417 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
418 cd write (iout,*) "Gather vbld_inv"
420 call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),
421 & MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),
422 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
423 cd write (iout,*) "Gather vbld side chain"
425 call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),
426 & MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),
427 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
428 cd write (iout,*) "Gather vbld_inv side chain"
430 call MPI_Allgatherv(vbld_inv(iint_start+nres),
431 & iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),
432 & iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
433 cd write (iout,*) "Gather theta"
435 call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),
436 & MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),
437 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
438 cd write (iout,*) "Gather phi"
440 call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),
441 & MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),
442 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
444 cd write (iout,*) "Gather alph"
446 call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),
447 & MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),
448 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
449 cd write (iout,*) "Gather omeg"
451 call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),
452 & MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),
453 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
455 time_gather=time_gather+MPI_Wtime()-time00
460 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
465 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
470 write (iout,1212) restyp(itype(i)),i,vbld(i),
471 &rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),
472 &rad2deg*alph(i),rad2deg*omeg(i)
475 1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
477 time_intfcart=time_intfcart+MPI_Wtime()-time01
481 c----------------------------------------------------------------------------
482 subroutine check_eint
483 C Check the gradient of energy in internal coordinates.
486 include 'COMMON.CONTROL'
487 include 'COMMON.CHAIN'
488 include 'COMMON.DERIV'
489 include 'COMMON.IOUNITS'
494 double precision x(maxvar),gana(maxvar),gg(maxvar)
496 double precision urparm(1)
497 double precision energia(0:n_ene),energia1(0:n_ene),
500 double precision fdum
502 double precision funcgrad,ff
505 double precision xi,etot,etot1,etot2
508 print '("Calling CHECK_INT",1pd12.3)',aincr
509 write (iout,'("Calling CHECK_INT",1pd12.3)') aincr
513 call geom_to_var(nvar,x)
514 call var_to_geom(nvar,x)
518 call etotal(energia(0))
520 call enerprint(energia(0))
523 if (MyID.ne.BossID) then
524 call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
532 cd write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
533 c write (iout,*) "Before gradient"
538 call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
540 c write (iout,*) "After gradient"
542 cd write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
547 call var_to_geom(nvar,x)
548 call chainbuild_extconf
549 call etotal(energia1(0))
552 call var_to_geom(nvar,x)
553 call chainbuild_extconf
554 call etotal(energia2(0))
556 gg(i)=(etot2-etot1)/aincr
557 c write (iout,*) i,etot1,etot2
560 write (iout,'(/2a)')' Variable Numerical Analytical',
566 else if (i.le.nphi+ntheta) then
569 else if (i.le.nphi+ntheta+nside) then
573 ii=i-(nphi+ntheta+nside)
576 write (iout,'(i3,a,i3,3(1pd16.6))')
577 & i,key,ii,gg(i),gana(i),
578 & 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)