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)
282 c call checkintcartgrad
285 write(iout,*) 'Calling CHECK_ECARTINT.'
288 call geom_to_var(nvar,x)
289 if (.not.split_ene) then
290 call etotal(energia(0))
292 call enerprint(energia(0))
294 write (iout,*) "enter cartgrad"
297 write (iout,*) "exit cartgrad"
301 write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
304 grad_s(j,0)=gcart(j,0)
308 grad_s(j,i)=gcart(j,i)
309 grad_s(j+3,i)=gxcart(j,i)
313 !- split gradient check
314 write (iout,*) "split_ene not supported"
316 c call etotal_long(energia(0))
317 c call enerprint(energia(0))
319 c write (iout,*) "enter cartgrad"
322 c write (iout,*) "exit cartgrad"
325 c write (iout,*) "longrange grad"
327 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
328 c & (gxcart(j,i),j=1,3)
331 c grad_s(j,0)=gcart(j,0)
335 c grad_s(j,i)=gcart(j,i)
336 c grad_s(j+3,i)=gxcart(j,i)
340 c call etotal_short(energia(0))
341 c call enerprint(energia(0))
343 c write (iout,*) "enter cartgrad"
346 c write (iout,*) "exit cartgrad"
349 c write (iout,*) "shortrange grad"
351 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
352 c & (gxcart(j,i),j=1,3)
355 c grad_s1(j,0)=gcart(j,0)
359 c grad_s1(j,i)=gcart(j,i)
360 c grad_s1(j+3,i)=gxcart(j,i)
364 write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
371 dcnorm_safe(k)=dc_norm(k,i)
372 dxnorm_safe(k)=dc_norm(k,i+nres)
379 c Broadcast the order to compute internal coordinates to the slaves.
381 c & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
383 c call int_from_cart1(.false.)
384 if (.not.split_ene) then
385 call etotal(energia1(0))
389 c call etotal_long(energia1(0))
391 c call etotal_short(energia1(0))
393 c write (iout,*) "etot11",etot11," etot12",etot12
395 !- end split gradient
396 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
399 c call int_from_cart1(.false.)
400 if (.not.split_ene) then
401 call etotal(energia1(0))
403 ggg(j)=(etot1-etot2)/(2*aincr)
406 c call etotal_long(energia1(0))
408 c ggg(j)=(etot11-etot21)/(2*aincr)
409 c call etotal_short(energia1(0))
411 c ggg1(j)=(etot12-etot22)/(2*aincr)
412 !- end split gradient
413 c write (iout,*) "etot21",etot21," etot22",etot22
415 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
420 dc(j,i+nres)=ddx(j)+aincr
422 c write (iout,*) "i",i," j",j," dxnorm+ and dxnorm"
423 c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
424 c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
425 c write (iout,*) "dxnormnorm",dsqrt(
426 c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
427 c write (iout,*) "dxnormnormsafe",dsqrt(
428 c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
430 if (.not.split_ene) then
431 call etotal(energia1(0))
435 c call etotal_long(energia1(0))
437 c call etotal_short(energia1(0))
440 !- end split gradient
441 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
442 dc(j,i+nres)=ddx(j)-aincr
444 c write (iout,*) "i",i," j",j," dxnorm- and dxnorm"
445 c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
446 c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
448 c write (iout,*) "dxnormnorm",dsqrt(
449 c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
450 c write (iout,*) "dxnormnormsafe",dsqrt(
451 c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
452 if (.not.split_ene) then
453 call etotal(energia1(0))
455 ggg(j+3)=(etot1-etot2)/(2*aincr)
458 c call etotal_long(energia1(0))
460 c ggg(j+3)=(etot11-etot21)/(2*aincr)
461 c call etotal_short(energia1(0))
463 c ggg1(j+3)=(etot12-etot22)/(2*aincr)
464 !- end split gradient
466 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
470 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
471 & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
473 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
474 & i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),
476 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
477 & i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),
478 & ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
483 c-------------------------------------------------------------------------
484 subroutine int_from_cart1(lprn)
485 implicit real*8 (a-h,o-z)
491 include 'COMMON.IOUNITS'
493 include 'COMMON.CHAIN'
495 include 'COMMON.INTERACT'
496 include 'COMMON.LOCAL'
497 include 'COMMON.NAMES'
498 include 'COMMON.SETUP'
499 include 'COMMON.TIME1'
501 if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
505 #if defined(PARINT) && defined(MPI)
506 do i=iint_start,iint_end
513 c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1
514 & +(c(j,i+1)-c(j,i))/dnorm2)
518 if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
519 if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
520 tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
522 if (itype(i-1).ne.10) then
523 tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
524 omicron(1,i)=alpha(i-2,i-1,i-1+nres)
525 omicron(2,i)=alpha(i-1+nres,i-1,i)
527 if (itype(i).ne.10) then
528 tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
531 omeg(i)=beta(nres+i,i,maxres2,i+1)
532 alph(i)=alpha(nres+i,i,maxres2)
533 theta(i+1)=alpha(i-1,i,i+1)
535 vbld_inv(i)=1.0d0/vbld(i)
536 vbld(nres+i)=dist(nres+i,i)
537 if (itype(i).ne.10) then
538 vbld_inv(nres+i)=1.0d0/vbld(nres+i)
540 vbld_inv(nres+i)=0.0d0
543 #if defined(PARINT) && defined(MPI)
544 if (nfgtasks1.gt.1) then
545 cd write(iout,*) "iint_start",iint_start," iint_count",
546 cd & (iint_count(i),i=0,nfgtasks-1)," iint_displ",
547 cd & (iint_displ(i),i=0,nfgtasks-1)
548 cd write (iout,*) "Gather vbld backbone"
551 call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),
552 & MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),
553 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
554 cd write (iout,*) "Gather vbld_inv"
556 call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),
557 & MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),
558 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
559 cd write (iout,*) "Gather vbld side chain"
561 call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),
562 & MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),
563 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
564 cd write (iout,*) "Gather vbld_inv side chain"
566 call MPI_Allgatherv(vbld_inv(iint_start+nres),
567 & iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),
568 & iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
569 cd write (iout,*) "Gather theta"
571 call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),
572 & MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),
573 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
574 cd write (iout,*) "Gather phi"
576 call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),
577 & MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),
578 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
580 cd write (iout,*) "Gather alph"
582 call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),
583 & MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),
584 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
585 cd write (iout,*) "Gather omeg"
587 call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),
588 & MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),
589 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
591 time_gather=time_gather+MPI_Wtime()-time00
596 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
601 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
606 write (iout,1212) restyp(itype(i)),i,vbld(i),
607 &rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),
608 &rad2deg*alph(i),rad2deg*omeg(i)
611 1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
613 time_intfcart=time_intfcart+MPI_Wtime()-time01
617 c----------------------------------------------------------------------------
618 subroutine check_eint
619 C Check the gradient of energy in internal coordinates.
620 implicit real*8 (a-h,o-z)
622 include 'COMMON.CHAIN'
623 include 'COMMON.DERIV'
624 include 'COMMON.IOUNITS'
628 dimension x(maxvar),gana(maxvar),gg(maxvar)
630 double precision urparm(1)
631 double precision energia(0:n_ene),energia1(0:n_ene),
637 print '(a)','Calling CHECK_INT.'
641 call geom_to_var(nvar,x)
642 call var_to_geom(nvar,x)
646 call etotal(energia(0))
648 call enerprint(energia(0))
651 if (MyID.ne.BossID) then
652 call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
660 cd write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
661 call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
662 cd write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
667 call var_to_geom(nvar,x)
669 call etotal(energia1(0))
672 call var_to_geom(nvar,x)
674 call etotal(energia2(0))
676 gg(i)=(etot2-etot1)/aincr
677 write (iout,*) i,etot1,etot2
680 write (iout,'(/2a)')' Variable Numerical Analytical',
686 else if (i.le.nphi+ntheta) then
689 else if (i.le.nphi+ntheta+nside) then
693 ii=i-(nphi+ntheta+nside)
696 write (iout,'(i3,a,i3,3(1pd16.6))')
697 & i,key,ii,gg(i),gana(i),
698 & 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)