1 C----------------------------------------------------------------------------
3 C Check the gradient of the energy in Cartesian coordinates.
6 include 'COMMON.CONTROL'
9 include 'COMMON.IOUNITS'
11 include 'COMMON.CONTACTS'
15 double precision ggg(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),
16 & g(maxvar),grad_s(6,maxres)
17 double precision energia(0:n_ene),energia1(0:n_ene)
18 double precision aincr2,etot,etot1,etot2
19 double precision dist,alpha,beta
20 double precision funcgrad,ff
24 double precision urparm(1)
31 print '("Calling CHECK_ECART",1pd12.3)',aincr
32 write (iout,'("Calling CHECK_ECART",1pd12.3)') aincr
35 call geom_to_var(nvar,x)
36 call etotal(energia(0))
38 call enerprint(energia(0))
42 call gradient(nvar,x,nf,g,uiparm,urparm,fdum)
46 write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3)
50 grad_s(j,i)=gradc(j,i,icg)
51 grad_s(j+3,i)=gradx(j,i,icg)
55 write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
66 c(j,k+nres)=c(j,k+nres)+aincr
68 call etotal(energia1(0))
70 ggg(j)=(etot1-etot)/aincr
74 c(j,k+nres)=c(j,k+nres)-aincr
78 c(j,i+nres)=c(j,i+nres)+aincr
79 dc(j,i+nres)=dc(j,i+nres)+aincr
80 call etotal(energia1(0))
82 ggg(j+3)=(etot1-etot)/aincr
86 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/)')
87 & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6)
91 c----------------------------------------------------------------------------
92 subroutine check_ecartint
93 C Check the gradient of the energy in Cartesian coordinates.
96 include 'COMMON.CONTROL'
97 include 'COMMON.CHAIN'
98 include 'COMMON.DERIV'
99 include 'COMMON.IOUNITS'
101 include 'COMMON.CONTACTS'
103 include 'COMMON.LOCAL'
104 include 'COMMON.SPLITELE'
107 double precision ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3),
108 & x(maxvar),g(maxvar)
109 double precision dcnorm_safe(3),dxnorm_safe(3)
110 double precision grad_s(6,0:maxres),grad_s1(6,0:maxres)
111 double precision phi_temp(maxres),theta_temp(maxres),
112 & alph_temp(maxres),omeg_temp(maxres)
113 double precision energia(0:n_ene),energia1(0:n_ene)
115 double precision urparm(1)
118 double precision etot,etot1,etot2,etot11,etot12,etot21,etot22
119 double precision dist,alpha,beta
126 call int_from_cart1(.false.)
129 c call checkintcartgrad
133 print '("Calling CHECK_ECARTINT",1pd12.3)',aincr
134 write (iout,'("Calling CHECK_ECARTINT",1pd12.3)') aincr
138 call geom_to_var(nvar,x)
139 if (.not.split_ene) then
140 call etotal(energia(0))
142 call enerprint(energia(0))
143 c write (iout,*) "enter cartgrad"
146 c write (iout,*) "exit cartgrad"
149 write (iout,'(//27(1h*)," Checking energy gradient ",27(1h*))')
150 write (iout,'(//4x,3a12,3x,3a12)')"gcart_x","gcart_y","gcart_z",
151 & "gxcart_x","gxcart_y","gxcart_z"
153 write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
154 & (gxcart(j,i),j=1,3)
157 grad_s(j,0)=gcart(j,0)
161 grad_s(j,i)=gcart(j,i)
162 grad_s(j+3,i)=gxcart(j,i)
166 !- split gradient check
168 call etotal_long(energia(0))
169 call enerprint(energia(0))
171 c write (iout,*) "enter cartgrad"
174 c write (iout,*) "exit cartgrad"
177 write (iout,*) "longrange grad"
179 write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
180 & (gxcart(j,i),j=1,3)
183 grad_s(j,0)=gcart(j,0)
187 grad_s(j,i)=gcart(j,i)
188 grad_s(j+3,i)=gxcart(j,i)
192 call etotal_short(energia(0))
193 call enerprint(energia(0))
195 c write (iout,*) "enter cartgrad"
198 c write (iout,*) "exit cartgrad"
201 write (iout,*) "shortrange grad"
203 write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3),
204 & (gxcart(j,i),j=1,3)
207 grad_s1(j,0)=gcart(j,0)
211 grad_s1(j,i)=gcart(j,i)
212 grad_s1(j+3,i)=gxcart(j,i)
216 c write (iout,*) "Vector dc"
218 c write (iout,'(i5,2(3f10.5,5x))')
219 c & i,(dc(j,i),j=1,3),(dc(j,i+nres),j=1,3)
221 c write (iout,*) "Coordinates after chainbuild_cart"
222 c call chainbuild_cart
224 write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
232 dcnorm_safe(k)=dc_norm(k,i)
233 dxnorm_safe(k)=dc_norm(k,i+nres)
240 c Broadcast the order to compute internal coordinates to the slaves.
242 c & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR)
244 c call int_from_cart1(.false.)
245 if (.not.split_ene) then
246 call etotal(energia1(0))
250 call etotal_long(energia1(0))
252 call etotal_short(energia1(0))
254 c write (iout,*) "etot11",etot11," etot12",etot12
256 !- end split gradient
257 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
261 c call int_from_cart1(.false.)
262 if (.not.split_ene) then
263 call etotal(energia1(0))
265 ggg(j)=(etot1-etot2)/(2*aincr)
268 call etotal_long(energia1(0))
270 ggg(j)=(etot11-etot21)/(2*aincr)
271 call etotal_short(energia1(0))
273 ggg1(j)=(etot12-etot22)/(2*aincr)
274 !- end split gradient
275 c write (iout,*) "etot21",etot21," etot22",etot22
277 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
282 dc(j,i+nres)=ddx(j)+aincr
284 c write (iout,*) "i",i," j",j," dxnorm+ and dxnorm"
285 c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
286 c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
287 c write (iout,*) "dxnormnorm",dsqrt(
288 c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
289 c write (iout,*) "dxnormnormsafe",dsqrt(
290 c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
292 if (.not.split_ene) then
293 call etotal(energia1(0))
297 call etotal_long(energia1(0))
299 call etotal_short(energia1(0))
302 !- end split gradient
303 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1
304 dc(j,i+nres)=ddx(j)-aincr
306 c write (iout,*) "i",i," j",j," dxnorm- and dxnorm"
307 c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3)
308 c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3)
310 c write (iout,*) "dxnormnorm",dsqrt(
311 c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2)
312 c write (iout,*) "dxnormnormsafe",dsqrt(
313 c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
314 if (.not.split_ene) then
315 call etotal(energia1(0))
317 ggg(j+3)=(etot1-etot2)/(2*aincr)
320 call etotal_long(energia1(0))
322 ggg(j+3)=(etot11-etot21)/(2*aincr)
323 call etotal_short(energia1(0))
325 ggg1(j+3)=(etot12-etot22)/(2*aincr)
326 !- end split gradient
328 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2
332 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
333 & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6)
335 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
336 & i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i),
338 write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)')
339 & i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6),
340 & ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6)
345 c-------------------------------------------------------------------------
346 subroutine int_from_cart1(lprn)
353 include 'COMMON.IOUNITS'
355 include 'COMMON.CHAIN'
357 include 'COMMON.INTERACT'
358 include 'COMMON.LOCAL'
359 include 'COMMON.NAMES'
360 include 'COMMON.SETUP'
361 include 'COMMON.TIME1'
364 double precision dnorm1,dnorm2,be
365 double precision time00
366 double precision dist,alpha,beta
367 double precision time01
368 if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
372 #if defined(PARINT) && defined(MPI)
373 do i=iint_start,iint_end
380 C print *,i,dnorm1,dnorm2
382 c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1
383 & +(c(j,i+1)-c(j,i))/dnorm2)
387 if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1)
388 if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then
389 tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
391 if (itype(i-1).ne.10) then
392 tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
393 omicron(1,i)=alpha(i-2,i-1,i-1+nres)
394 omicron(2,i)=alpha(i-1+nres,i-1,i)
396 if (itype(i).ne.10) then
397 tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
400 omeg(i)=beta(nres+i,i,maxres2,i+1)
402 alph(i)=alpha(nres+i,i,maxres2)
404 theta(i+1)=alpha(i-1,i,i+1)
407 vbld_inv(i)=1.0d0/vbld(i)
408 vbld(nres+i)=dist(nres+i,i)
409 C print *,vbld(i+nres)
411 if (itype(i).ne.10) then
412 vbld_inv(nres+i)=1.0d0/vbld(nres+i)
414 vbld_inv(nres+i)=0.0d0
417 #if defined(PARINT) && defined(MPI)
418 if (nfgtasks1.gt.1) then
419 cd write(iout,*) "iint_start",iint_start," iint_count",
420 cd & (iint_count(i),i=0,nfgtasks-1)," iint_displ",
421 cd & (iint_displ(i),i=0,nfgtasks-1)
422 cd write (iout,*) "Gather vbld backbone"
425 call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1),
426 & MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0),
427 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
428 cd write (iout,*) "Gather vbld_inv"
430 call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1),
431 & MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0),
432 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
433 cd write (iout,*) "Gather vbld side chain"
435 call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1),
436 & MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0),
437 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
438 cd write (iout,*) "Gather vbld_inv side chain"
440 call MPI_Allgatherv(vbld_inv(iint_start+nres),
441 & iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1),
442 & iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
443 cd write (iout,*) "Gather theta"
445 call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1),
446 & MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0),
447 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
448 cd write (iout,*) "Gather phi"
450 call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1),
451 & MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0),
452 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
454 cd write (iout,*) "Gather alph"
456 call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1),
457 & MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0),
458 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
459 cd write (iout,*) "Gather omeg"
461 call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1),
462 & MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0),
463 & MPI_DOUBLE_PRECISION,FG_COMM1,IERR)
465 time_gather=time_gather+MPI_Wtime()-time00
470 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
475 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
480 write (iout,1212) restyp(itype(i)),i,vbld(i),
481 &rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),
482 &rad2deg*alph(i),rad2deg*omeg(i)
485 1212 format (a3,'(',i3,')',2(f15.10,2f10.2))
487 time_intfcart=time_intfcart+MPI_Wtime()-time01
491 c----------------------------------------------------------------------------
492 subroutine check_eint
493 C Check the gradient of energy in internal coordinates.
496 include 'COMMON.CONTROL'
497 include 'COMMON.CHAIN'
498 include 'COMMON.DERIV'
499 include 'COMMON.IOUNITS'
504 double precision x(maxvar),gana(maxvar),gg(maxvar)
506 double precision urparm(1)
507 double precision energia(0:n_ene),energia1(0:n_ene),
510 double precision fdum
512 double precision funcgrad,ff
515 double precision xi,etot,etot1,etot2
518 print '("Calling CHECK_INT",1pd12.3)',aincr
519 write (iout,'("Calling CHECK_INT",1pd12.3)') aincr
523 call geom_to_var(nvar,x)
524 call var_to_geom(nvar,x)
528 call etotal(energia(0))
530 call enerprint(energia(0))
533 if (MyID.ne.BossID) then
534 call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID)
542 cd write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar)
543 c write (iout,*) "Before gradient"
548 call gradient(nvar,x,nf,gana,uiparm,urparm,fdum)
550 c write (iout,*) "After gradient"
552 cd write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar)
557 call var_to_geom(nvar,x)
558 call chainbuild_extconf
559 call etotal(energia1(0))
562 call var_to_geom(nvar,x)
563 call chainbuild_extconf
564 call etotal(energia2(0))
566 gg(i)=(etot2-etot1)/aincr
567 c write (iout,*) i,etot1,etot2
570 write (iout,'(/2a)')' Variable Numerical Analytical',
576 else if (i.le.nphi+ntheta) then
579 else if (i.le.nphi+ntheta+nside) then
583 ii=i-(nphi+ntheta+nside)
586 write (iout,'(i3,a,i3,3(1pd16.6))')
587 & i,key,ii,gg(i),gana(i),
588 & 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr)