1 subroutine check_cartgrad
2 C Check the gradient of Cartesian coordinates in internal coordinates.
3 implicit real*8 (a-h,o-z)
5 include 'COMMON.IOUNITS'
10 include 'COMMON.DERIV'
11 include 'COMMON.SCCOR'
12 dimension temp(6,maxres),xx(3),gg(3)
13 indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1
15 * Check the gradient of the virtual-bond and SC vectors in the internal
21 write (iout,'(a)') '**************** dx/dalpha'
27 temp(k,i)=dc(k,nres+i)
31 gg(k)=(dc(k,nres+i)-temp(k,i))/aincr
32 xx(k)=dabs((gg(k)-dxds(k,i))/(aincr*dabs(dxds(k,i))+aincr))
34 write (iout,'(i4,3e15.6/4x,3e15.6,3f9.3)')
35 & i,(gg(k),k=1,3),(dxds(k,i),k=1,3),(xx(k),k=1,3)
41 write (iout,'(a)') '**************** dx/domega'
47 temp(k,i)=dc(k,nres+i)
51 gg(k)=(dc(k,nres+i)-temp(k,i))/aincr
52 xx(k)=dabs((gg(k)-dxds(k+3,i))/
53 & (aincr*dabs(dxds(k+3,i))+aincr))
55 write (iout,'(i4,3e15.6/4x,3e15.6,3f9.3)')
56 & i,(gg(k),k=1,3),(dxds(k+3,i),k=1,3),(xx(k),k=1,3)
62 write (iout,'(a)') '**************** dx/dtheta'
66 theta(i)=theta(i)+aincr
69 temp(k,j)=dc(k,nres+j)
75 c print *,'i=',i-2,' j=',j-1,' ii=',ii
77 gg(k)=(dc(k,nres+j)-temp(k,j))/aincr
78 xx(k)=dabs((gg(k)-dxdv(k,ii))/
79 & (aincr*dabs(dxdv(k,ii))+aincr))
81 write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)')
82 & i,j,(gg(k),k=1,3),(dxdv(k,ii),k=1,3),(xx(k),k=1,3)
89 write (iout,'(a)') '***************** dx/dphi'
95 temp(k,j)=dc(k,nres+j)
103 gg(k)=(dc(k,nres+j)-temp(k,j))/aincr
104 xx(k)=dabs((gg(k)-dxdv(k+3,ii))/
105 & (aincr*dabs(dxdv(k+3,ii))+aincr))
107 write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)')
108 & i,j,(gg(k),k=1,3),(dxdv(k+3,ii),k=1,3),(xx(k),k=1,3)
114 write (iout,'(a)') '****************** ddc/dtheta'
117 theta(i+2)=thet+aincr
128 gg(k)=(dc(k,j)-temp(k,j))/aincr
129 xx(k)=dabs((gg(k)-dcdv(k,ii))/
130 & (aincr*dabs(dcdv(k,ii))+aincr))
132 write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)')
133 & i,j,(gg(k),k=1,3),(dcdv(k,ii),k=1,3),(xx(k),k=1,3)
143 write (iout,'(a)') '******************* ddc/dphi'
157 gg(k)=(dc(k,j)-temp(k,j))/aincr
158 xx(k)=dabs((gg(k)-dcdv(k+3,ii))/
159 & (aincr*dabs(dcdv(k+3,ii))+aincr))
161 write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)')
162 & i,j,(gg(k),k=1,3),(dcdv(k+3,ii),k=1,3),(xx(k),k=1,3)
174 C----------------------------------------------------------------------------
175 subroutine check_ecart
176 C Check the gradient of the energy in Cartesian coordinates.
177 implicit real*8 (a-h,o-z)
179 include 'COMMON.CHAIN'
180 include 'COMMON.DERIV'
181 include 'COMMON.IOUNITS'
183 include 'COMMON.CONTACTS'
184 include 'COMMON.SCCOR'
186 dimension ggg(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),g(maxvar)
187 dimension grad_s(6,maxres)
188 double precision energia(0:n_ene),energia1(0:n_ene)
190 double precision urparm(1)
197 print '(a)','CG processor',me,' calling CHECK_CART.'
200 call geom_to_var(nvar,x)
201 call etotal(energia(0))
203 call enerprint(energia(0))
204 call gradient(nvar,x,nf,g,uiparm,urparm,fdum)
207 write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
211 grad_s(j,i)=gradc(j,i,icg)
212 grad_s(j+3,i)=gradx(j,i,icg)
216 write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
224 dc(j,i)=dc(j,i)+aincr
227 c(j,k+nres)=c(j,k+nres)+aincr
229 call etotal(energia1(0))
231 ggg(j)=(etot1-etot)/aincr
235 c(j,k+nres)=c(j,k+nres)-aincr
239 c(j,i+nres)=c(j,i+nres)+aincr
240 dc(j,i+nres)=dc(j,i+nres)+aincr
241 call etotal(energia1(0))
243 ggg(j+3)=(etot1-etot)/aincr
247 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/)')
248 & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6)
252 c----------------------------------------------------------------------------
253 subroutine check_ecartint
254 C Check the gradient of the energy in Cartesian coordinates.
255 implicit real*8 (a-h,o-z)
257 include 'COMMON.CONTROL'
258 include 'COMMON.CHAIN'
259 include 'COMMON.DERIV'
260 include 'COMMON.IOUNITS'
262 include 'COMMON.CONTACTS'
264 include 'COMMON.LOCAL'
265 include 'COMMON.SPLITELE'
266 include 'COMMON.SCCOR'
268 dimension ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),
270 dimension dcnorm_safe(3),dxnorm_safe(3)
271 dimension grad_s(6,0:maxres),grad_s1(6,0:maxres)
272 double precision phi_temp(maxres),theta_temp(maxres),
273 & alph_temp(maxres),omeg_temp(maxres)
274 double precision energia(0:n_ene),energia1(0:n_ene)
276 double precision urparm(1)
285 c call checkintcartgrad
288 write(iout,*) 'Calling CHECK_ECARTINT.'
291 call geom_to_var(nvar,x)
292 if (.not.split_ene) then
293 call etotal(energia(0))
295 c write (iout,*) "atu?", gloc_sc(1,i,icg),gloc(i,icg)
298 call enerprint(energia(0))
300 write (iout,*) "enter cartgrad"
302 c write (iout,*) gloc_sc(1,i,icg)
306 write (iout,*) "exit cartgrad"
310 write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
313 grad_s(j,0)=gcart(j,0)
317 grad_s(j,i)=gcart(j,i)
318 grad_s(j+3,i)=gxcart(j,i)
322 !- split gradient check
324 call etotal_long(energia(0))
325 call enerprint(energia(0))
327 write (iout,*) "enter cartgrad"
330 write (iout,*) "exit cartgrad"
333 write (iout,*) "longrange grad"
335 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
336 & (gxcart(j,i),j=1,3)
339 grad_s(j,0)=gcart(j,0)
343 grad_s(j,i)=gcart(j,i)
344 grad_s(j+3,i)=gxcart(j,i)
348 call etotal_short(energia(0))
349 call enerprint(energia(0))
351 c write (iout,*) gloc_sc(1,i,icg)
354 write (iout,*) "enter cartgrad"
357 write (iout,*) "exit cartgrad"
360 write (iout,*) "shortrange grad"
362 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
363 & (gxcart(j,i),j=1,3)
366 grad_s1(j,0)=gcart(j,0)
370 grad_s1(j,i)=gcart(j,i)
371 grad_s1(j+3,i)=gxcart(j,i)
375 write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
382 dcnorm_safe(k)=dc_norm(k,i)
383 dxnorm_safe(k)=dc_norm(k,i+nres)
390 c Broadcast the order to compute internal coordinates to the slaves.
392 c & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
394 c call int_from_cart1(.false.)
395 if (.not.split_ene) then
396 call etotal(energia1(0))
400 call etotal_long(energia1(0))
402 call etotal_short(energia1(0))
404 c write (iout,*) "etot11",etot11," etot12",etot12
406 !- end split gradient
407 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
410 c call int_from_cart1(.false.)
411 if (.not.split_ene) then
412 call etotal(energia1(0))
414 ggg(j)=(etot1-etot2)/(2*aincr)
417 call etotal_long(energia1(0))
419 ggg(j)=(etot11-etot21)/(2*aincr)
420 call etotal_short(energia1(0))
422 ggg1(j)=(etot12-etot22)/(2*aincr)
423 !- end split gradient
424 c write (iout,*) "etot21",etot21," etot22",etot22
426 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
431 dc(j,i+nres)=ddx(j)+aincr
433 c write (iout,*) "i",i," j",j," dxnorm+ and dxnorm"
434 c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
435 c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
436 c write (iout,*) "dxnormnorm",dsqrt(
437 c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
438 c write (iout,*) "dxnormnormsafe",dsqrt(
439 c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
441 if (.not.split_ene) then
442 call etotal(energia1(0))
446 call etotal_long(energia1(0))
448 call etotal_short(energia1(0))
451 !- end split gradient
452 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
453 dc(j,i+nres)=ddx(j)-aincr
455 c write (iout,*) "i",i," j",j," dxnorm- and dxnorm"
456 c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
457 c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
459 c write (iout,*) "dxnormnorm",dsqrt(
460 c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
461 c write (iout,*) "dxnormnormsafe",dsqrt(
462 c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
463 if (.not.split_ene) then
464 call etotal(energia1(0))
466 ggg(j+3)=(etot1-etot2)/(2*aincr)
469 call etotal_long(energia1(0))
471 ggg(j+3)=(etot11-etot21)/(2*aincr)
472 call etotal_short(energia1(0))
474 ggg1(j+3)=(etot12-etot22)/(2*aincr)
475 !- end split gradient
477 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
481 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
482 & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
484 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
485 & i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),
487 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
488 & i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),
489 & ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
494 c-------------------------------------------------------------------------
495 subroutine int_from_cart1(lprn)
496 implicit real*8 (a-h,o-z)
502 include 'COMMON.IOUNITS'
504 include 'COMMON.CHAIN'
506 include 'COMMON.INTERACT'
507 include 'COMMON.LOCAL'
508 include 'COMMON.NAMES'
509 include 'COMMON.SETUP'
510 include 'COMMON.TIME1'
512 if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
516 #if defined(PARINT) && defined(MPI)
517 do i=iint_start,iint_end+1
524 c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1
525 & +(c(j,i+1)-c(j,i))/dnorm2)
529 if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
530 if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
531 tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
533 if (itype(i-1).ne.10) then
534 tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
535 omicron(1,i)=alpha(i-2,i-1,i-1+nres)
536 omicron(2,i)=alpha(i-1+nres,i-1,i)
538 if (itype(i).ne.10) then
539 tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
542 omeg(i)=beta(nres+i,i,maxres2,i+1)
543 alph(i)=alpha(nres+i,i,maxres2)
544 theta(i+1)=alpha(i-1,i,i+1)
546 vbld_inv(i)=1.0d0/vbld(i)
547 vbld(nres+i)=dist(nres+i,i)
548 if (itype(i).ne.10) then
549 vbld_inv(nres+i)=1.0d0/vbld(nres+i)
551 vbld_inv(nres+i)=0.0d0
555 #if defined(PARINT) && defined(MPI)
556 if (nfgtasks1.gt.1) then
557 cd write(iout,*) "iint_start",iint_start," iint_count",
558 cd & (iint_count(i),i=0,nfgtasks-1)," iint_displ",
559 cd & (iint_displ(i),i=0,nfgtasks-1)
560 cd write (iout,*) "Gather vbld backbone"
563 call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),
564 & MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),
565 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
566 cd write (iout,*) "Gather vbld_inv"
568 call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),
569 & MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),
570 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
571 cd write (iout,*) "Gather vbld side chain"
573 call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),
574 & MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),
575 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
576 cd write (iout,*) "Gather vbld_inv side chain"
578 call MPI_Allgatherv(vbld_inv(iint_start+nres),
579 & iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),
580 & iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
581 cd write (iout,*) "Gather theta"
583 call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),
584 & MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),
585 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
586 cd write (iout,*) "Gather phi"
588 call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),
589 & MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),
590 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
592 cd write (iout,*) "Gather alph"
594 call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),
595 & MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),
596 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
597 cd write (iout,*) "Gather omeg"
599 call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),
600 & MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),
601 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
603 time_gather=time_gather+MPI_Wtime()-time00
608 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
613 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
618 write (iout,1212) restyp(itype(i)),i,vbld(i),
619 &rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),
620 &rad2deg*alph(i),rad2deg*omeg(i)
623 1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
625 time_intfcart=time_intfcart+MPI_Wtime()-time01
629 c----------------------------------------------------------------------------
630 subroutine check_eint
631 C Check the gradient of energy in internal coordinates.
632 implicit real*8 (a-h,o-z)
634 include 'COMMON.CHAIN'
635 include 'COMMON.DERIV'
636 include 'COMMON.IOUNITS'
640 dimension x(maxvar),gana(maxvar),gg(maxvar)
642 double precision urparm(1)
643 double precision energia(0:n_ene),energia1(0:n_ene),
649 print '(a)','Calling CHECK_INT.'
653 call geom_to_var(nvar,x)
654 call var_to_geom(nvar,x)
658 call etotal(energia(0))
660 call enerprint(energia(0))
663 if (MyID.ne.BossID) then
664 call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
672 cd write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
673 call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
674 cd write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
679 call var_to_geom(nvar,x)
681 call etotal(energia1(0))
684 call var_to_geom(nvar,x)
686 call etotal(energia2(0))
688 gg(i)=(etot2-etot1)/aincr
689 write (iout,*) i,etot1,etot2
692 write (iout,'(/2a)')' Variable Numerical Analytical',
698 else if (i.le.nphi+ntheta) then
701 else if (i.le.nphi+ntheta+nside) then
705 ii=i-(nphi+ntheta+nside)
708 write (iout,'(i3,a,i3,3(1pd16.6))')
709 & i,key,ii,gg(i),gana(i),
710 & 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)