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
315 call etotal_long(energia(0))
316 call enerprint(energia(0))
318 write (iout,*) "enter cartgrad"
321 write (iout,*) "exit cartgrad"
324 write (iout,*) "longrange grad"
326 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
327 & (gxcart(j,i),j=1,3)
330 grad_s(j,0)=gcart(j,0)
334 grad_s(j,i)=gcart(j,i)
335 grad_s(j+3,i)=gxcart(j,i)
339 call etotal_short(energia(0))
340 call enerprint(energia(0))
342 write (iout,*) "enter cartgrad"
345 write (iout,*) "exit cartgrad"
348 write (iout,*) "shortrange grad"
350 write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
351 & (gxcart(j,i),j=1,3)
354 grad_s1(j,0)=gcart(j,0)
358 grad_s1(j,i)=gcart(j,i)
359 grad_s1(j+3,i)=gxcart(j,i)
363 write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
370 dcnorm_safe(k)=dc_norm(k,i)
371 dxnorm_safe(k)=dc_norm(k,i+nres)
378 c Broadcast the order to compute internal coordinates to the slaves.
380 c & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
382 c call int_from_cart1(.false.)
383 if (.not.split_ene) then
384 call etotal(energia1(0))
388 call etotal_long(energia1(0))
390 call etotal_short(energia1(0))
392 c write (iout,*) "etot11",etot11," etot12",etot12
394 !- end split gradient
395 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
398 c call int_from_cart1(.false.)
399 if (.not.split_ene) then
400 call etotal(energia1(0))
402 ggg(j)=(etot1-etot2)/(2*aincr)
405 call etotal_long(energia1(0))
407 ggg(j)=(etot11-etot21)/(2*aincr)
408 call etotal_short(energia1(0))
410 ggg1(j)=(etot12-etot22)/(2*aincr)
411 !- end split gradient
412 c write (iout,*) "etot21",etot21," etot22",etot22
414 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
419 dc(j,i+nres)=ddx(j)+aincr
421 c write (iout,*) "i",i," j",j," dxnorm+ and dxnorm"
422 c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
423 c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
424 c write (iout,*) "dxnormnorm",dsqrt(
425 c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
426 c write (iout,*) "dxnormnormsafe",dsqrt(
427 c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
429 if (.not.split_ene) then
430 call etotal(energia1(0))
434 call etotal_long(energia1(0))
436 call etotal_short(energia1(0))
439 !- end split gradient
440 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
441 dc(j,i+nres)=ddx(j)-aincr
443 c write (iout,*) "i",i," j",j," dxnorm- and dxnorm"
444 c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
445 c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
447 c write (iout,*) "dxnormnorm",dsqrt(
448 c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
449 c write (iout,*) "dxnormnormsafe",dsqrt(
450 c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
451 if (.not.split_ene) then
452 call etotal(energia1(0))
454 ggg(j+3)=(etot1-etot2)/(2*aincr)
457 call etotal_long(energia1(0))
459 ggg(j+3)=(etot11-etot21)/(2*aincr)
460 call etotal_short(energia1(0))
462 ggg1(j+3)=(etot12-etot22)/(2*aincr)
463 !- end split gradient
465 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
469 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
470 & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
472 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
473 & i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),
475 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
476 & i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),
477 & ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
482 c-------------------------------------------------------------------------
483 subroutine int_from_cart1(lprn)
484 implicit real*8 (a-h,o-z)
490 include 'COMMON.IOUNITS'
492 include 'COMMON.CHAIN'
494 include 'COMMON.INTERACT'
495 include 'COMMON.LOCAL'
496 include 'COMMON.NAMES'
497 include 'COMMON.SETUP'
498 include 'COMMON.TIME1'
500 if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
504 #if defined(PARINT) && defined(MPI)
505 do i=iint_start,iint_end
512 c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1
513 & +(c(j,i+1)-c(j,i))/dnorm2)
516 if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
517 omeg(i)=beta(nres+i,i,maxres2,i+1)
518 alph(i)=alpha(nres+i,i,maxres2)
519 theta(i+1)=alpha(i-1,i,i+1)
521 vbld_inv(i)=1.0d0/vbld(i)
522 vbld(nres+i)=dist(nres+i,i)
523 if (itype(i).ne.10) then
524 vbld_inv(nres+i)=1.0d0/vbld(nres+i)
526 vbld_inv(nres+i)=0.0d0
529 #if defined(PARINT) && defined(MPI)
530 if (nfgtasks1.gt.1) then
531 cd write(iout,*) "iint_start",iint_start," iint_count",
532 cd & (iint_count(i),i=0,nfgtasks-1)," iint_displ",
533 cd & (iint_displ(i),i=0,nfgtasks-1)
534 cd write (iout,*) "Gather vbld backbone"
537 call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),
538 & MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),
539 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
540 cd write (iout,*) "Gather vbld_inv"
542 call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),
543 & MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),
544 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
545 cd write (iout,*) "Gather vbld side chain"
547 call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),
548 & MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),
549 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
550 cd write (iout,*) "Gather vbld_inv side chain"
552 call MPI_Allgatherv(vbld_inv(iint_start+nres),
553 & iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),
554 & iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
555 cd write (iout,*) "Gather theta"
557 call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),
558 & MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),
559 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
560 cd write (iout,*) "Gather phi"
562 call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),
563 & MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),
564 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
566 cd write (iout,*) "Gather alph"
568 call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),
569 & MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),
570 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
571 cd write (iout,*) "Gather omeg"
573 call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),
574 & MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),
575 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
577 time_gather=time_gather+MPI_Wtime()-time00
582 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
587 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
592 write (iout,1212) restyp(itype(i)),i,vbld(i),
593 &rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),
594 &rad2deg*alph(i),rad2deg*omeg(i)
597 1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
599 time_intfcart=time_intfcart+MPI_Wtime()-time01
603 c----------------------------------------------------------------------------
604 subroutine check_eint
605 C Check the gradient of energy in internal coordinates.
606 implicit real*8 (a-h,o-z)
608 include 'COMMON.CHAIN'
609 include 'COMMON.DERIV'
610 include 'COMMON.IOUNITS'
614 dimension x(maxvar),gana(maxvar),gg(maxvar)
616 double precision urparm(1)
617 double precision energia(0:n_ene),energia1(0:n_ene),
623 print '(a)','Calling CHECK_INT.'
627 call geom_to_var(nvar,x)
628 call var_to_geom(nvar,x)
632 call etotal(energia(0))
634 call enerprint(energia(0))
637 if (MyID.ne.BossID) then
638 call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
646 cd write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
647 call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
648 cd write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
653 call var_to_geom(nvar,x)
655 call etotal(energia1(0))
658 call var_to_geom(nvar,x)
660 call etotal(energia2(0))
662 gg(i)=(etot2-etot1)/aincr
663 write (iout,*) i,etot1,etot2
666 write (iout,'(/2a)')' Variable Numerical Analytical',
672 else if (i.le.nphi+ntheta) then
675 else if (i.le.nphi+ntheta+nside) then
679 ii=i-(nphi+ntheta+nside)
682 write (iout,'(i3,a,i3,3(1pd16.6))')
683 & i,key,ii,gg(i),gana(i),
684 & 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)