C---------------------------------------------------------------------------- subroutine check_ecart C Check the gradient of the energy in Cartesian coordinates. implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.IOUNITS' include 'COMMON.VAR' include 'COMMON.CONTACTS' integer i,j,k integer icall common /srutu/ icall double precision ggg(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar), & g(maxvar),grad_s(6,maxres) double precision energia(0:n_ene),energia1(0:n_ene) double precision aincr2,etot,etot1,etot2 double precision dist,alpha,beta double precision funcgrad,ff external funcgrad integer nf integer uiparm(1) double precision urparm(1) double precision fdum external fdum icg=1 nf=0 nfl=0 call zerograd print '("Calling CHECK_ECART",1pd12.3)',aincr write (iout,'("Calling CHECK_ECART",1pd12.3)') aincr nf=0 icall=0 call geom_to_var(nvar,x) call etotal(energia(0)) etot=energia(0) call enerprint(energia(0)) #ifdef LBFGS ff=funcgrad(x,g) #else call gradient(nvar,x,nf,g,uiparm,urparm,fdum) #endif icall =1 do i=1,nres write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3) enddo do i=1,nres do j=1,3 grad_s(j,i)=gradc(j,i,icg) grad_s(j+3,i)=gradx(j,i,icg) enddo enddo call flush(iout) write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors' do i=1,nres do j=1,3 xx(j)=c(j,i+nres) ddc(j)=dc(j,i) ddx(j)=dc(j,i+nres) enddo do j=1,3 dc(j,i)=dc(j,i)+aincr do k=i+1,nres c(j,k)=c(j,k)+aincr c(j,k+nres)=c(j,k+nres)+aincr enddo call etotal(energia1(0)) etot1=energia1(0) ggg(j)=(etot1-etot)/aincr dc(j,i)=ddc(j) do k=i+1,nres c(j,k)=c(j,k)-aincr c(j,k+nres)=c(j,k+nres)-aincr enddo enddo do j=1,3 c(j,i+nres)=c(j,i+nres)+aincr dc(j,i+nres)=dc(j,i+nres)+aincr call etotal(energia1(0)) etot1=energia1(0) ggg(j+3)=(etot1-etot)/aincr c(j,i+nres)=xx(j) dc(j,i+nres)=ddx(j) enddo write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/)') & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6) enddo return end c---------------------------------------------------------------------------- subroutine check_ecartint C Check the gradient of the energy in Cartesian coordinates. implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.IOUNITS' include 'COMMON.VAR' include 'COMMON.CONTACTS' include 'COMMON.MD' include 'COMMON.LOCAL' include 'COMMON.SPLITELE' integer icall common /srutu/ icall double precision ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3), & x(maxvar),g(maxvar) double precision dcnorm_safe(3),dxnorm_safe(3) double precision grad_s(6,0:maxres),grad_s1(6,0:maxres) double precision phi_temp(maxres),theta_temp(maxres), & alph_temp(maxres),omeg_temp(maxres) double precision energia(0:n_ene),energia1(0:n_ene) integer uiparm(1) double precision urparm(1) external fdum integer i,j,k,nf double precision etot,etot1,etot2,etot11,etot12,etot21,etot22 double precision dist,alpha,beta c r_cut=2.0d0 c rlambd=0.3d0 icg=1 nf=0 nfl=0 c print *,"ATU 3" call int_from_cart1(.false.) call intout c call intcartderiv c call checkintcartgrad call zerograd c aincr=8.0D-7 c aincr=1.0D-7 print '("Calling CHECK_ECARTINT",1pd12.3)',aincr write (iout,'("Calling CHECK_ECARTINT",1pd12.3)') aincr call cartprint nf=0 icall=0 call geom_to_var(nvar,x) if (.not.split_ene) then call etotal(energia(0)) etot=energia(0) call enerprint(energia(0)) c write (iout,*) "enter cartgrad" c call flush(iout) call cartgrad c write (iout,*) "exit cartgrad" c call flush(iout) icall =1 write (iout,'(//27(1h*)," Checking energy gradient ",27(1h*))') write (iout,'(//4x,3a12,3x,3a12)')"gcart_x","gcart_y","gcart_z", & "gxcart_x","gxcart_y","gxcart_z" do i=1,nres write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3), & (gxcart(j,i),j=1,3) enddo do j=1,3 grad_s(j,0)=gcart(j,0) enddo do i=1,nres do j=1,3 grad_s(j,i)=gcart(j,i) grad_s(j+3,i)=gxcart(j,i) enddo enddo else !- split gradient check call zerograd call etotal_long(energia(0)) call enerprint(energia(0)) call flush(iout) c write (iout,*) "enter cartgrad" c call flush(iout) call cartgrad c write (iout,*) "exit cartgrad" c call flush(iout) icall =1 write (iout,*) "longrange grad" do i=1,nres write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3), & (gxcart(j,i),j=1,3) enddo do j=1,3 grad_s(j,0)=gcart(j,0) enddo do i=1,nres do j=1,3 grad_s(j,i)=gcart(j,i) grad_s(j+3,i)=gxcart(j,i) enddo enddo call zerograd call etotal_short(energia(0)) call enerprint(energia(0)) call flush(iout) c write (iout,*) "enter cartgrad" c call flush(iout) call cartgrad c write (iout,*) "exit cartgrad" c call flush(iout) icall =1 write (iout,*) "shortrange grad" do i=1,nres write (iout,'(i4,3e12.4,3x,3e12.4)') i,(gcart(j,i),j=1,3), & (gxcart(j,i),j=1,3) enddo do j=1,3 grad_s1(j,0)=gcart(j,0) enddo do i=1,nres do j=1,3 grad_s1(j,i)=gcart(j,i) grad_s1(j+3,i)=gxcart(j,i) enddo enddo endif c write (iout,*) "Vector dc" c do i=0,nres c write (iout,'(i5,2(3f10.5,5x))') c & i,(dc(j,i),j=1,3),(dc(j,i+nres),j=1,3) c enddo c write (iout,*) "Coordinates after chainbuild_cart" c call chainbuild_cart c call cartprint write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors' do i=0,nres print *,i do j=1,3 xx(j)=c(j,i+nres) ddc(j)=dc(j,i) ddx(j)=dc(j,i+nres) do k=1,3 dcnorm_safe(k)=dc_norm(k,i) dxnorm_safe(k)=dc_norm(k,i+nres) enddo enddo do j=1,3 dc(j,i)=ddc(j)+aincr call chainbuild_cart #ifdef MPI c Broadcast the order to compute internal coordinates to the slaves. c if (nfgtasks.gt.1) c & call MPI_Bcast(6,1,MPI_INTEGER,king,FG_COMM,IERROR) #endif c call int_from_cart1(.false.) if (.not.split_ene) then call etotal(energia1(0)) etot1=energia1(0) else !- split gradient call etotal_long(energia1(0)) etot11=energia1(0) call etotal_short(energia1(0)) etot12=energia1(0) c write (iout,*) "etot11",etot11," etot12",etot12 endif !- end split gradient c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1 dc(j,i)=ddc(j)-aincr call chainbuild_cart C print *,c(j,i) c call int_from_cart1(.false.) if (.not.split_ene) then call etotal(energia1(0)) etot2=energia1(0) ggg(j)=(etot1-etot2)/(2*aincr) else !- split gradient call etotal_long(energia1(0)) etot21=energia1(0) ggg(j)=(etot11-etot21)/(2*aincr) call etotal_short(energia1(0)) etot22=energia1(0) ggg1(j)=(etot12-etot22)/(2*aincr) !- end split gradient c write (iout,*) "etot21",etot21," etot22",etot22 endif c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2 dc(j,i)=ddc(j) call chainbuild_cart enddo do j=1,3 dc(j,i+nres)=ddx(j)+aincr call chainbuild_cart c write (iout,*) "i",i," j",j," dxnorm+ and dxnorm" c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3) c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3) c write (iout,*) "dxnormnorm",dsqrt( c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2) c write (iout,*) "dxnormnormsafe",dsqrt( c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2) c write (iout,*) if (.not.split_ene) then call etotal(energia1(0)) etot1=energia1(0) else !- split gradient call etotal_long(energia1(0)) etot11=energia1(0) call etotal_short(energia1(0)) etot12=energia1(0) endif !- end split gradient c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1 dc(j,i+nres)=ddx(j)-aincr call chainbuild_cart c write (iout,*) "i",i," j",j," dxnorm- and dxnorm" c write (iout,'(3f15.10)') (dc_norm(k,i+nres),k=1,3) c write (iout,'(3f15.10)') (dxnorm_safe(k),k=1,3) c write (iout,*) c write (iout,*) "dxnormnorm",dsqrt( c & dc_norm(1,i+nres)**2+dc_norm(2,i+nres)**2+dc_norm(3,i+nres)**2) c write (iout,*) "dxnormnormsafe",dsqrt( c & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2) if (.not.split_ene) then call etotal(energia1(0)) etot2=energia1(0) ggg(j+3)=(etot1-etot2)/(2*aincr) else !- split gradient call etotal_long(energia1(0)) etot21=energia1(0) ggg(j+3)=(etot11-etot21)/(2*aincr) call etotal_short(energia1(0)) etot22=energia1(0) ggg1(j+3)=(etot12-etot22)/(2*aincr) !- end split gradient endif c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2 dc(j,i+nres)=ddx(j) call chainbuild_cart enddo write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') & i,(ggg(k),k=1,6),(grad_s(k,i),k=1,6),(ggg(k)/grad_s(k,i),k=1,6) if (split_ene) then write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') & i,(ggg1(k),k=1,6),(grad_s1(k,i),k=1,6),(ggg1(k)/grad_s1(k,i), & k=1,6) write (iout,'(i3,6(1pe12.5)/3x,6(1pe12.5)/3x,6(1pe12.5)/)') & i,(ggg(k)+ggg1(k),k=1,6),(grad_s(k,i)+grad_s1(k,i),k=1,6), & ((ggg(k)+ggg1(k))/(grad_s(k,i)+grad_s1(k,i)),k=1,6) endif enddo return end c------------------------------------------------------------------------- subroutine int_from_cart1(lprn) implicit none include 'DIMENSIONS' #ifdef MPI include 'mpif.h' integer ierror #endif include 'COMMON.IOUNITS' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.GEO' include 'COMMON.INTERACT' include 'COMMON.LOCAL' include 'COMMON.NAMES' include 'COMMON.SETUP' include 'COMMON.TIME1' logical lprn integer i,j double precision dnorm1,dnorm2,be double precision time00 double precision dist,alpha,beta double precision time01 if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates' #ifdef TIMING time01=MPI_Wtime() #endif #if defined(PARINT) && defined(MPI) do i=iint_start,iint_end #else do i=2,nres #endif C print *,i dnorm1=dist(i-1,i) dnorm2=dist(i,i+1) C print *,i,dnorm1,dnorm2 do j=1,3 c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1 & +(c(j,i+1)-c(j,i))/dnorm2) enddo be=0.0D0 if (i.gt.2) then if (i.le.nres) phi(i+1)=beta(i-2,i-1,i,i+1) if ((itype(i).ne.10).and.(itype(i-1).ne.10)) then tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres) endif if (itype(i-1).ne.10) then tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1) omicron(1,i)=alpha(i-2,i-1,i-1+nres) omicron(2,i)=alpha(i-1+nres,i-1,i) endif if (itype(i).ne.10) then tauangle(2,i+1)=beta(i-2,i-1,i,i+nres) endif endif omeg(i)=beta(nres+i,i,maxres2,i+1) C print *,omeg(i) alph(i)=alpha(nres+i,i,maxres2) C print *,alph(i) theta(i+1)=alpha(i-1,i,i+1) vbld(i)=dist(i-1,i) C print *,vbld(i) vbld_inv(i)=1.0d0/vbld(i) vbld(nres+i)=dist(nres+i,i) C print *,vbld(i+nres) if (itype(i).ne.10) then vbld_inv(nres+i)=1.0d0/vbld(nres+i) else vbld_inv(nres+i)=0.0d0 endif enddo #if defined(PARINT) && defined(MPI) if (nfgtasks1.gt.1) then cd write(iout,*) "iint_start",iint_start," iint_count", cd & (iint_count(i),i=0,nfgtasks-1)," iint_displ", cd & (iint_displ(i),i=0,nfgtasks-1) cd write (iout,*) "Gather vbld backbone" cd call flush(iout) time00=MPI_Wtime() call MPI_Allgatherv(vbld(iint_start),iint_count(fg_rank1), & MPI_DOUBLE_PRECISION,vbld(1),iint_count(0),iint_displ(0), & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) cd write (iout,*) "Gather vbld_inv" cd call flush(iout) call MPI_Allgatherv(vbld_inv(iint_start),iint_count(fg_rank1), & MPI_DOUBLE_PRECISION,vbld_inv(1),iint_count(0),iint_displ(0), & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) cd write (iout,*) "Gather vbld side chain" cd call flush(iout) call MPI_Allgatherv(vbld(iint_start+nres),iint_count(fg_rank1), & MPI_DOUBLE_PRECISION,vbld(nres+1),iint_count(0),iint_displ(0), & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) cd write (iout,*) "Gather vbld_inv side chain" cd call flush(iout) call MPI_Allgatherv(vbld_inv(iint_start+nres), & iint_count(fg_rank1),MPI_DOUBLE_PRECISION,vbld_inv(nres+1), & iint_count(0),iint_displ(0),MPI_DOUBLE_PRECISION,FG_COMM1,IERR) cd write (iout,*) "Gather theta" cd call flush(iout) call MPI_Allgatherv(theta(iint_start+1),iint_count(fg_rank1), & MPI_DOUBLE_PRECISION,theta(2),iint_count(0),iint_displ(0), & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) cd write (iout,*) "Gather phi" cd call flush(iout) call MPI_Allgatherv(phi(iint_start+1),iint_count(fg_rank1), & MPI_DOUBLE_PRECISION,phi(2),iint_count(0),iint_displ(0), & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) #ifdef CRYST_SC cd write (iout,*) "Gather alph" cd call flush(iout) call MPI_Allgatherv(alph(iint_start),iint_count(fg_rank1), & MPI_DOUBLE_PRECISION,alph(1),iint_count(0),iint_displ(0), & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) cd write (iout,*) "Gather omeg" cd call flush(iout) call MPI_Allgatherv(omeg(iint_start),iint_count(fg_rank1), & MPI_DOUBLE_PRECISION,omeg(1),iint_count(0),iint_displ(0), & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) #endif time_gather=time_gather+MPI_Wtime()-time00 endif #endif do i=1,nres-1 do j=1,3 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1) enddo enddo do i=2,nres-1 do j=1,3 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres) enddo enddo if (lprn) then do i=2,nres write (iout,1212) restyp(itype(i)),i,vbld(i), &rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i), &rad2deg*alph(i),rad2deg*omeg(i) enddo endif 1212 format (a3,'(',i3,')',2(f15.10,2f10.2)) #ifdef TIMING time_intfcart=time_intfcart+MPI_Wtime()-time01 #endif return end c---------------------------------------------------------------------------- subroutine check_eint C Check the gradient of energy in internal coordinates. implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.IOUNITS' include 'COMMON.VAR' include 'COMMON.GEO' integer icall common /srutu/ icall double precision x(maxvar),gana(maxvar),gg(maxvar) integer uiparm(1) double precision urparm(1) double precision energia(0:n_ene),energia1(0:n_ene), & energia2(0:n_ene) character*6 key double precision fdum external fdum double precision funcgrad,ff external funcgrad integer i,ii,nf double precision xi,etot,etot1,etot2 call zerograd c aincr=1.0D-7 print '("Calling CHECK_INT",1pd12.3)',aincr write (iout,'("Calling CHECK_INT",1pd12.3)') aincr nf=0 nfl=0 icg=1 call geom_to_var(nvar,x) call var_to_geom(nvar,x) call chainbuild icall=1 print *,'ICG=',ICG call etotal(energia(0)) etot = energia(0) call enerprint(energia(0)) print *,'ICG=',ICG #ifdef MPL if (MyID.ne.BossID) then call mp_bcast(x(1),8*(nvar+3),BossID,fgGroupID) nf=x(nvar+1) nfl=x(nvar+2) icg=x(nvar+3) endif #endif nf=1 nfl=3 cd write (iout,'(10f8.3)') (rad2deg*x(i),i=1,nvar) c write (iout,*) "Before gradient" c call flush(iout) #ifdef LBFGS ff=funcgrad(x,gana) #else call gradient(nvar,x,nf,gana,uiparm,urparm,fdum) #endif c write (iout,*) "After gradient" c call flush(iout) cd write (iout,'(i3,1pe14.4)') (i,gana(i),i=1,nvar) icall=1 do i=1,nvar xi=x(i) x(i)=xi-0.5D0*aincr call var_to_geom(nvar,x) call chainbuild_extconf call etotal(energia1(0)) etot1=energia1(0) x(i)=xi+0.5D0*aincr call var_to_geom(nvar,x) call chainbuild_extconf call etotal(energia2(0)) etot2=energia2(0) gg(i)=(etot2-etot1)/aincr c write (iout,*) i,etot1,etot2 x(i)=xi enddo write (iout,'(/2a)')' Variable Numerical Analytical', & ' RelDiff*100% ' do i=1,nvar if (i.le.nphi) then ii=i key = ' phi' else if (i.le.nphi+ntheta) then ii=i-nphi key=' theta' else if (i.le.nphi+ntheta+nside) then ii=i-(nphi+ntheta) key=' alpha' else ii=i-(nphi+ntheta+nside) key=' omega' endif write (iout,'(i3,a,i3,3(1pd16.6))') & i,key,ii,gg(i),gana(i), & 100.0D0*dabs(gg(i)-gana(i))/(dabs(gana(i))+aincr) enddo return end