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
400 c call int_from_cart1(.false.)
401 if (.not.split_ene) then
402 call etotal(energia1(0))
404 ggg(j)=(etot1-etot2)/(2*aincr)
407 call etotal_long(energia1(0))
409 ggg(j)=(etot11-etot21)/(2*aincr)
410 call etotal_short(energia1(0))
412 ggg1(j)=(etot12-etot22)/(2*aincr)
413 !- end split gradient
414 c write (iout,*) "etot21",etot21," etot22",etot22
416 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
421 dc(j,i+nres)=ddx(j)+aincr
423 c write (iout,*) "i",i," j",j," dxnorm+ and dxnorm"
424 c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
425 c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
426 c write (iout,*) "dxnormnorm",dsqrt(
427 c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
428 c write (iout,*) "dxnormnormsafe",dsqrt(
429 c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
431 if (.not.split_ene) then
432 call etotal(energia1(0))
436 call etotal_long(energia1(0))
438 call etotal_short(energia1(0))
441 !- end split gradient
442 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
443 dc(j,i+nres)=ddx(j)-aincr
445 c write (iout,*) "i",i," j",j," dxnorm- and dxnorm"
446 c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
447 c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
449 c write (iout,*) "dxnormnorm",dsqrt(
450 c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
451 c write (iout,*) "dxnormnormsafe",dsqrt(
452 c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
453 if (.not.split_ene) then
454 call etotal(energia1(0))
456 ggg(j+3)=(etot1-etot2)/(2*aincr)
459 call etotal_long(energia1(0))
461 ggg(j+3)=(etot11-etot21)/(2*aincr)
462 call etotal_short(energia1(0))
464 ggg1(j+3)=(etot12-etot22)/(2*aincr)
465 !- end split gradient
467 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
471 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
472 & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
474 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
475 & i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),
477 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
478 & i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),
479 & ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
484 c-------------------------------------------------------------------------
485 subroutine int_from_cart1(lprn)
486 implicit real*8 (a-h,o-z)
492 include 'COMMON.IOUNITS'
494 include 'COMMON.CHAIN'
496 include 'COMMON.INTERACT'
497 include 'COMMON.LOCAL'
498 include 'COMMON.NAMES'
499 include 'COMMON.SETUP'
500 include 'COMMON.TIME1'
502 if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
506 #if defined(PARINT) && defined(MPI)
507 do i=iint_start,iint_end
514 C print *,i,dnorm1,dnorm2
516 c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1
517 & +(c(j,i+1)-c(j,i))/dnorm2)
521 if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
522 if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
523 tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
525 if (itype(i-1).ne.10) then
526 tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
527 omicron(1,i)=alpha(i-2,i-1,i-1+nres)
528 omicron(2,i)=alpha(i-1+nres,i-1,i)
530 if (itype(i).ne.10) then
531 tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
534 omeg(i)=beta(nres+i,i,maxres2,i+1)
536 alph(i)=alpha(nres+i,i,maxres2)
538 theta(i+1)=alpha(i-1,i,i+1)
541 vbld_inv(i)=1.0d0/vbld(i)
542 vbld(nres+i)=dist(nres+i,i)
543 C print *,vbld(i+nres)
545 if (itype(i).ne.10) then
546 vbld_inv(nres+i)=1.0d0/vbld(nres+i)
548 vbld_inv(nres+i)=0.0d0
551 #if defined(PARINT) && defined(MPI)
552 if (nfgtasks1.gt.1) then
553 cd write(iout,*) "iint_start",iint_start," iint_count",
554 cd & (iint_count(i),i=0,nfgtasks-1)," iint_displ",
555 cd & (iint_displ(i),i=0,nfgtasks-1)
556 cd write (iout,*) "Gather vbld backbone"
559 call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),
560 & MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),
561 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
562 cd write (iout,*) "Gather vbld_inv"
564 call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),
565 & MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),
566 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
567 cd write (iout,*) "Gather vbld side chain"
569 call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),
570 & MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),
571 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
572 cd write (iout,*) "Gather vbld_inv side chain"
574 call MPI_Allgatherv(vbld_inv(iint_start+nres),
575 & iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),
576 & iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
577 cd write (iout,*) "Gather theta"
579 call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),
580 & MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),
581 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
582 cd write (iout,*) "Gather phi"
584 call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),
585 & MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),
586 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
588 cd write (iout,*) "Gather alph"
590 call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),
591 & MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),
592 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
593 cd write (iout,*) "Gather omeg"
595 call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),
596 & MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),
597 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
599 time_gather=time_gather+MPI_Wtime()-time00
604 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
609 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
614 write (iout,1212) restyp(itype(i)),i,vbld(i),
615 &rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),
616 &rad2deg*alph(i),rad2deg*omeg(i)
619 1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
621 time_intfcart=time_intfcart+MPI_Wtime()-time01
625 c----------------------------------------------------------------------------
626 subroutine check_eint
627 C Check the gradient of energy in internal coordinates.
628 implicit real*8 (a-h,o-z)
630 include 'COMMON.CHAIN'
631 include 'COMMON.DERIV'
632 include 'COMMON.IOUNITS'
636 dimension x(maxvar),gana(maxvar),gg(maxvar)
638 double precision urparm(1)
639 double precision energia(0:n_ene),energia1(0:n_ene),
645 print '(a)','Calling CHECK_INT.'
649 call geom_to_var(nvar,x)
650 call var_to_geom(nvar,x)
654 call etotal(energia(0))
656 call enerprint(energia(0))
659 if (MyID.ne.BossID) then
660 call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
668 cd write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
669 call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
670 cd write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
675 call var_to_geom(nvar,x)
677 call etotal(energia1(0))
680 call var_to_geom(nvar,x)
682 call etotal(energia2(0))
684 gg(i)=(etot2-etot1)/aincr
685 write (iout,*) i,etot1,etot2
688 write (iout,'(/2a)')' Variable Numerical Analytical',
694 else if (i.le.nphi+ntheta) then
697 else if (i.le.nphi+ntheta+nside) then
701 ii=i-(nphi+ntheta+nside)
704 write (iout,'(i3,a,i3,3(1pd16.6))')
705 & i,key,ii,gg(i),gana(i),
706 & 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)