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 dimension temp(6,maxres),xx(3),gg(3)
12 indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1
14 * Check the gradient of the virtual-bond and SC vectors in the internal
20 write (iout,'(a)') '**************** dx/dalpha'
26 temp(k,i)=dc(k,nres+i)
30 gg(k)=(dc(k,nres+i)-temp(k,i))/aincr
31 xx(k)=dabs((gg(k)-dxds(k,i))/(aincr*dabs(dxds(k,i))+aincr))
33 write (iout,'(i4,3e15.6/4x,3e15.6,3f9.3)')
34 & i,(gg(k),k=1,3),(dxds(k,i),k=1,3),(xx(k),k=1,3)
40 write (iout,'(a)') '**************** dx/domega'
46 temp(k,i)=dc(k,nres+i)
50 gg(k)=(dc(k,nres+i)-temp(k,i))/aincr
51 xx(k)=dabs((gg(k)-dxds(k+3,i))/
52 & (aincr*dabs(dxds(k+3,i))+aincr))
54 write (iout,'(i4,3e15.6/4x,3e15.6,3f9.3)')
55 & i,(gg(k),k=1,3),(dxds(k+3,i),k=1,3),(xx(k),k=1,3)
61 write (iout,'(a)') '**************** dx/dtheta'
65 theta(i)=theta(i)+aincr
68 temp(k,j)=dc(k,nres+j)
74 c print *,'i=',i-2,' j=',j-1,' ii=',ii
76 gg(k)=(dc(k,nres+j)-temp(k,j))/aincr
77 xx(k)=dabs((gg(k)-dxdv(k,ii))/
78 & (aincr*dabs(dxdv(k,ii))+aincr))
80 write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)')
81 & i,j,(gg(k),k=1,3),(dxdv(k,ii),k=1,3),(xx(k),k=1,3)
88 write (iout,'(a)') '***************** dx/dphi'
94 temp(k,j)=dc(k,nres+j)
102 gg(k)=(dc(k,nres+j)-temp(k,j))/aincr
103 xx(k)=dabs((gg(k)-dxdv(k+3,ii))/
104 & (aincr*dabs(dxdv(k+3,ii))+aincr))
106 write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)')
107 & i,j,(gg(k),k=1,3),(dxdv(k+3,ii),k=1,3),(xx(k),k=1,3)
113 write (iout,'(a)') '****************** ddc/dtheta'
116 theta(i+2)=thet+aincr
127 gg(k)=(dc(k,j)-temp(k,j))/aincr
128 xx(k)=dabs((gg(k)-dcdv(k,ii))/
129 & (aincr*dabs(dcdv(k,ii))+aincr))
131 write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)')
132 & i,j,(gg(k),k=1,3),(dcdv(k,ii),k=1,3),(xx(k),k=1,3)
142 write (iout,'(a)') '******************* ddc/dphi'
156 gg(k)=(dc(k,j)-temp(k,j))/aincr
157 xx(k)=dabs((gg(k)-dcdv(k+3,ii))/
158 & (aincr*dabs(dcdv(k+3,ii))+aincr))
160 write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)')
161 & i,j,(gg(k),k=1,3),(dcdv(k+3,ii),k=1,3),(xx(k),k=1,3)
173 C----------------------------------------------------------------------------
174 subroutine check_ecart
175 C Check the gradient of the energy in Cartesian coordinates.
176 implicit real*8 (a-h,o-z)
178 include 'COMMON.CHAIN'
179 include 'COMMON.DERIV'
180 include 'COMMON.IOUNITS'
182 include 'COMMON.CONTACTS'
184 dimension ggg(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),g(maxvar)
185 dimension grad_s(6,maxres)
186 double precision energia(0:n_ene),energia1(0:n_ene)
188 double precision urparm(1)
195 print '(a)','CG processor',me,' calling CHECK_CART.'
198 call geom_to_var(nvar,x)
199 call etotal(energia(0))
201 call enerprint(energia(0))
202 call gradient(nvar,x,nf,g,uiparm,urparm,fdum)
205 write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
209 grad_s(j,i)=gradc(j,i,icg)
210 grad_s(j+3,i)=gradx(j,i,icg)
214 write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
222 dc(j,i)=dc(j,i)+aincr
225 c(j,k+nres)=c(j,k+nres)+aincr
227 call etotal(energia1(0))
229 ggg(j)=(etot1-etot)/aincr
233 c(j,k+nres)=c(j,k+nres)-aincr
237 c(j,i+nres)=c(j,i+nres)+aincr
238 dc(j,i+nres)=dc(j,i+nres)+aincr
239 call etotal(energia1(0))
241 ggg(j+3)=(etot1-etot)/aincr
245 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/)')
246 & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6)
250 c----------------------------------------------------------------------------
251 subroutine check_ecartint
252 C Check the gradient of the energy in Cartesian coordinates.
253 implicit real*8 (a-h,o-z)
255 include 'COMMON.CONTROL'
256 include 'COMMON.CHAIN'
257 include 'COMMON.DERIV'
258 include 'COMMON.IOUNITS'
260 include 'COMMON.CONTACTS'
262 include 'COMMON.LOCAL'
263 include 'COMMON.SPLITELE'
265 dimension ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),
267 dimension dcnorm_safe(3),dxnorm_safe(3)
268 dimension grad_s(6,0:maxres),grad_s1(6,0:maxres)
269 double precision phi_temp(maxres),theta_temp(maxres),
270 & alph_temp(maxres),omeg_temp(maxres)
271 double precision energia(0:n_ene),energia1(0:n_ene)
273 double precision urparm(1)
283 c call checkintcartgrad
286 write(iout,*) 'Calling CHECK_ECARTINT.'
289 call geom_to_var(nvar,x)
290 if (.not.split_ene) then
291 call etotal(energia(0))
293 call enerprint(energia(0))
295 write (iout,*) "enter cartgrad"
298 write (iout,*) "exit cartgrad"
302 write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
305 grad_s(j,0)=gcart(j,0)
309 grad_s(j,i)=gcart(j,i)
310 grad_s(j+3,i)=gxcart(j,i)
314 !- split gradient check
316 call etotal_long(energia(0))
317 call enerprint(energia(0))
319 write (iout,*) "enter cartgrad"
322 write (iout,*) "exit cartgrad"
325 write (iout,*) "longrange grad"
327 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
328 & (gxcart(j,i),j=1,3)
331 grad_s(j,0)=gcart(j,0)
335 grad_s(j,i)=gcart(j,i)
336 grad_s(j+3,i)=gxcart(j,i)
340 call etotal_short(energia(0))
341 call enerprint(energia(0))
343 write (iout,*) "enter cartgrad"
346 write (iout,*) "exit cartgrad"
349 write (iout,*) "shortrange grad"
351 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
352 & (gxcart(j,i),j=1,3)
355 grad_s1(j,0)=gcart(j,0)
359 grad_s1(j,i)=gcart(j,i)
360 grad_s1(j+3,i)=gxcart(j,i)
364 write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
372 dcnorm_safe(k)=dc_norm(k,i)
373 dxnorm_safe(k)=dc_norm(k,i+nres)
380 c Broadcast the order to compute internal coordinates to the slaves.
382 c & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
384 c call int_from_cart1(.false.)
385 if (.not.split_ene) then
386 call etotal(energia1(0))
390 call etotal_long(energia1(0))
392 call etotal_short(energia1(0))
394 c write (iout,*) "etot11",etot11," etot12",etot12
396 !- end split gradient
397 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
401 c call int_from_cart1(.false.)
402 if (.not.split_ene) then
403 call etotal(energia1(0))
405 ggg(j)=(etot1-etot2)/(2*aincr)
408 call etotal_long(energia1(0))
410 ggg(j)=(etot11-etot21)/(2*aincr)
411 call etotal_short(energia1(0))
413 ggg1(j)=(etot12-etot22)/(2*aincr)
414 !- end split gradient
415 c write (iout,*) "etot21",etot21," etot22",etot22
417 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
422 dc(j,i+nres)=ddx(j)+aincr
424 c write (iout,*) "i",i," j",j," dxnorm+ and dxnorm"
425 c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
426 c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
427 c write (iout,*) "dxnormnorm",dsqrt(
428 c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
429 c write (iout,*) "dxnormnormsafe",dsqrt(
430 c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
432 if (.not.split_ene) then
433 call etotal(energia1(0))
437 call etotal_long(energia1(0))
439 call etotal_short(energia1(0))
442 !- end split gradient
443 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
444 dc(j,i+nres)=ddx(j)-aincr
446 c write (iout,*) "i",i," j",j," dxnorm- and dxnorm"
447 c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
448 c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
450 c write (iout,*) "dxnormnorm",dsqrt(
451 c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
452 c write (iout,*) "dxnormnormsafe",dsqrt(
453 c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
454 if (.not.split_ene) then
455 call etotal(energia1(0))
457 ggg(j+3)=(etot1-etot2)/(2*aincr)
460 call etotal_long(energia1(0))
462 ggg(j+3)=(etot11-etot21)/(2*aincr)
463 call etotal_short(energia1(0))
465 ggg1(j+3)=(etot12-etot22)/(2*aincr)
466 !- end split gradient
468 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
472 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
473 & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
475 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
476 & i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),
478 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
479 & i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),
480 & ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
485 c-------------------------------------------------------------------------
486 subroutine int_from_cart1(lprn)
487 implicit real*8 (a-h,o-z)
493 include 'COMMON.IOUNITS'
495 include 'COMMON.CHAIN'
497 include 'COMMON.INTERACT'
498 include 'COMMON.LOCAL'
499 include 'COMMON.NAMES'
500 include 'COMMON.SETUP'
501 include 'COMMON.TIME1'
503 if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
507 #if defined(PARINT) && defined(MPI)
508 do i=iint_start,iint_end
515 C print *,i,dnorm1,dnorm2
517 c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1
518 & +(c(j,i+1)-c(j,i))/dnorm2)
522 if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
523 if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
524 tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
526 if (itype(i-1).ne.10) then
527 tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
528 omicron(1,i)=alpha(i-2,i-1,i-1+nres)
529 omicron(2,i)=alpha(i-1+nres,i-1,i)
531 if (itype(i).ne.10) then
532 tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
535 omeg(i)=beta(nres+i,i,maxres2,i+1)
537 alph(i)=alpha(nres+i,i,maxres2)
539 theta(i+1)=alpha(i-1,i,i+1)
542 vbld_inv(i)=1.0d0/vbld(i)
543 vbld(nres+i)=dist(nres+i,i)
544 C print *,vbld(i+nres)
546 if (itype(i).ne.10) then
547 vbld_inv(nres+i)=1.0d0/vbld(nres+i)
549 vbld_inv(nres+i)=0.0d0
552 #if defined(PARINT) && defined(MPI)
553 if (nfgtasks1.gt.1) then
554 cd write(iout,*) "iint_start",iint_start," iint_count",
555 cd & (iint_count(i),i=0,nfgtasks-1)," iint_displ",
556 cd & (iint_displ(i),i=0,nfgtasks-1)
557 cd write (iout,*) "Gather vbld backbone"
560 call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),
561 & MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),
562 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
563 cd write (iout,*) "Gather vbld_inv"
565 call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),
566 & MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),
567 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
568 cd write (iout,*) "Gather vbld side chain"
570 call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),
571 & MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),
572 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
573 cd write (iout,*) "Gather vbld_inv side chain"
575 call MPI_Allgatherv(vbld_inv(iint_start+nres),
576 & iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),
577 & iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
578 cd write (iout,*) "Gather theta"
580 call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),
581 & MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),
582 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
583 cd write (iout,*) "Gather phi"
585 call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),
586 & MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),
587 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
589 cd write (iout,*) "Gather alph"
591 call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),
592 & MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),
593 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
594 cd write (iout,*) "Gather omeg"
596 call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),
597 & MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),
598 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
600 time_gather=time_gather+MPI_Wtime()-time00
605 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
610 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
615 write (iout,1212) restyp(itype(i)),i,vbld(i),
616 &rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),
617 &rad2deg*alph(i),rad2deg*omeg(i)
620 1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
622 time_intfcart=time_intfcart+MPI_Wtime()-time01
626 c----------------------------------------------------------------------------
627 subroutine check_eint
628 C Check the gradient of energy in internal coordinates.
629 implicit real*8 (a-h,o-z)
631 include 'COMMON.CHAIN'
632 include 'COMMON.DERIV'
633 include 'COMMON.IOUNITS'
637 dimension x(maxvar),gana(maxvar),gg(maxvar)
639 double precision urparm(1)
640 double precision energia(0:n_ene),energia1(0:n_ene),
646 print '(a)','Calling CHECK_INT.'
650 call geom_to_var(nvar,x)
651 call var_to_geom(nvar,x)
655 call etotal(energia(0))
657 call enerprint(energia(0))
660 if (MyID.ne.BossID) then
661 call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
669 cd write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
670 call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
671 cd write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
676 call var_to_geom(nvar,x)
678 call etotal(energia1(0))
681 call var_to_geom(nvar,x)
683 call etotal(energia2(0))
685 gg(i)=(etot2-etot1)/aincr
686 write (iout,*) i,etot1,etot2
689 write (iout,'(/2a)')' Variable Numerical Analytical',
695 else if (i.le.nphi+ntheta) then
698 else if (i.le.nphi+ntheta+nside) then
702 ii=i-(nphi+ntheta+nside)
705 write (iout,'(i3,a,i3,3(1pd16.6))')
706 & i,key,ii,gg(i),gana(i),
707 & 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)