subroutine check_cartgrad C Check the gradient of Cartesian coordinates in internal coordinates. implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.IOUNITS' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.DERIV' double precision temp(6,maxres),xx(3),gg(3),thet,theti,phii,alphi, & omegi,aincr2 integer indmat integer i,ii,j,k indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1 integer nf * * Check the gradient of the virtual-bond and SC vectors in the internal * coordinates. * print '("Calling CHECK_ECART",1pd12.3)',aincr write (iout,'("Calling CHECK_ECART",1pd12.3)') aincr aincr2=0.5d0*aincr call chainbuild_extconf call cartder write (iout,'(a)') '**************** dx/dalpha' write (iout,'(a)') do i=2,nres-1 alphi=alph(i) alph(i)=alph(i)+aincr do k=1,3 temp(k,i)=dc(k,nres+i) enddo call chainbuild_extconf do k=1,3 gg(k)=(dc(k,nres+i)-temp(k,i))/aincr xx(k)=dabs((gg(k)-dxds(k,i))/(aincr*dabs(dxds(k,i))+aincr)) enddo write (iout,'(i4,3e15.6/4x,3e15.6,3f9.3)') & i,(gg(k),k=1,3),(dxds(k,i),k=1,3),(xx(k),k=1,3) write (iout,'(a)') alph(i)=alphi call chainbuild_extconf enddo write (iout,'(a)') write (iout,'(a)') '**************** dx/domega' write (iout,'(a)') do i=2,nres-1 omegi=omeg(i) omeg(i)=omeg(i)+aincr do k=1,3 temp(k,i)=dc(k,nres+i) enddo call chainbuild_extconf do k=1,3 gg(k)=(dc(k,nres+i)-temp(k,i))/aincr xx(k)=dabs((gg(k)-dxds(k+3,i))/ & (aincr*dabs(dxds(k+3,i))+aincr)) enddo write (iout,'(i4,3e15.6/4x,3e15.6,3f9.3)') & i,(gg(k),k=1,3),(dxds(k+3,i),k=1,3),(xx(k),k=1,3) write (iout,'(a)') omeg(i)=omegi call chainbuild_extconf enddo write (iout,'(a)') write (iout,'(a)') '**************** dx/dtheta' write (iout,'(a)') do i=3,nres theti=theta(i) theta(i)=theta(i)+aincr do j=i-1,nres-1 do k=1,3 temp(k,j)=dc(k,nres+j) enddo enddo call chainbuild_extconf do j=i-1,nres-1 ii = indmat(i-2,j) c print *,'i=',i-2,' j=',j-1,' ii=',ii do k=1,3 gg(k)=(dc(k,nres+j)-temp(k,j))/aincr xx(k)=dabs((gg(k)-dxdv(k,ii))/ & (aincr*dabs(dxdv(k,ii))+aincr)) enddo write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') & i,j,(gg(k),k=1,3),(dxdv(k,ii),k=1,3),(xx(k),k=1,3) write(iout,'(a)') enddo write (iout,'(a)') theta(i)=theti call chainbuild_extconf enddo write (iout,'(a)') '***************** dx/dphi' write (iout,'(a)') do i=4,nres phi(i)=phi(i)+aincr do j=i-1,nres-1 do k=1,3 temp(k,j)=dc(k,nres+j) enddo enddo call chainbuild_extconf do j=i-1,nres-1 ii = indmat(i-2,j) c print *,'ii=',ii do k=1,3 gg(k)=(dc(k,nres+j)-temp(k,j))/aincr xx(k)=dabs((gg(k)-dxdv(k+3,ii))/ & (aincr*dabs(dxdv(k+3,ii))+aincr)) enddo write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') & i,j,(gg(k),k=1,3),(dxdv(k+3,ii),k=1,3),(xx(k),k=1,3) write(iout,'(a)') enddo phi(i)=phi(i)-aincr call chainbuild_extconf enddo write (iout,'(a)') '****************** ddc/dtheta' do i=1,nres-2 thet=theta(i+2) theta(i+2)=thet+aincr do j=i,nres do k=1,3 temp(k,j)=dc(k,j) enddo enddo call chainbuild_extconf do j=i+1,nres-1 ii = indmat(i,j) c print *,'ii=',ii do k=1,3 gg(k)=(dc(k,j)-temp(k,j))/aincr xx(k)=dabs((gg(k)-dcdv(k,ii))/ & (aincr*dabs(dcdv(k,ii))+aincr)) enddo write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') & i,j,(gg(k),k=1,3),(dcdv(k,ii),k=1,3),(xx(k),k=1,3) write (iout,'(a)') enddo do j=1,nres do k=1,3 dc(k,j)=temp(k,j) enddo enddo theta(i+2)=thet enddo write (iout,'(a)') '******************* ddc/dphi' do i=1,nres-3 phii=phi(i+3) phi(i+3)=phii+aincr do j=1,nres do k=1,3 temp(k,j)=dc(k,j) enddo enddo call chainbuild_extconf do j=i+2,nres-1 ii = indmat(i+1,j) c print *,'ii=',ii do k=1,3 gg(k)=(dc(k,j)-temp(k,j))/aincr xx(k)=dabs((gg(k)-dcdv(k+3,ii))/ & (aincr*dabs(dcdv(k+3,ii))+aincr)) enddo write (iout,'(2i4,3e14.6/8x,3e14.6,3f9.3)') & i,j,(gg(k),k=1,3),(dcdv(k+3,ii),k=1,3),(xx(k),k=1,3) write (iout,'(a)') enddo do j=1,nres do k=1,3 dc(k,j)=temp(k,j) enddo enddo phi(i+3)=phii enddo return end 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 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)) call gradient(nvar,x,nf,g,uiparm,urparm,fdum) 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 #ifdef FIVEDIAG !----------------------------------------------------------------------------- subroutine check_ecartint ! Check the gradient of the energy in Cartesian coordinates. implicit none include 'DIMENSIONS' include 'COMMON.CONTROL' include 'COMMON.CHAIN' include 'COMMON.INTERACT' 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 ddc1(3),ddcn(3),dcnorm_safe1(3),dcnorm_safe2(3) double precision energia(0:n_ene),energia1(0:n_ene) integer uiparm(1) double precision urparm(1) double precision fdum external fdum integer i,j,k,nf double precision etot,etot1,etot2,etot11,etot12,etot21,etot22 double precision dist,alpha,beta icg=1 nf=0 nfl=0 call intout ! call intcartderiv ! call checkintcartgrad call zerograd aincr=1.0D-5 write(iout,*) 'Calling CHECK_ECARTINT.' nf=0 icall=0 write (iout,*) "Before geom_to_var" call geom_to_var(nvar,x) write (iout,*) "after geom_to_var" write (iout,*) "split_ene ",split_ene call flush(iout) if (.not.split_ene) then write(iout,*) 'Calling CHECK_ECARTINT if' call etotal(energia) !elwrite(iout,*) 'Calling CHECK_ECARTINT if' etot=energia(0) write (iout,*) "etot",etot call enerprint(energia(0)) call flush(iout) !el call enerprint(energia) !elwrite(iout,*) 'Calling CHECK_ECARTINT if' call flush(iout) write (iout,*) "enter cartgrad" call flush(iout) call cartgrad !elwrite(iout,*) 'Calling CHECK_ECARTINT if' write (iout,*) "exit cartgrad" call flush(iout) icall =1 do i=1,nres write (iout,'(i5,3f10.5)') i,(gradxorr(j,i),j=1,3) enddo do j=1,3 grad_s(j,0)=gcart(j,0) enddo !elwrite(iout,*) 'Calling CHECK_ECARTINT if' 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 write(iout,*) 'Calling CHECK_ECARTIN else.' !- split gradient check call zerograd call etotal_long(energia) call enerprint(energia(0)) !el call enerprint(energia) call flush(iout) write (iout,*) "enter cartgrad" call flush(iout) call cartgrad write (iout,*) "exit cartgrad" call flush(iout) icall =1 write (iout,*) "longrange grad" do i=1,nres write (iout,'(i5,3f10.5,5x,3f10.5)') 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) call enerprint(energia(0)) call flush(iout) write (iout,*) "enter cartgrad" call flush(iout) call cartgrad write (iout,*) "exit cartgrad" call flush(iout) icall =1 write (iout,*) "shortrange grad" do i=1,nres write (iout,'(i5,3f10.5,5x,3f10.5)') 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 write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors' ! do i=1,nres c do i=nnt,nct do i=1,nres do j=1,3 if (nnt.gt.1 .and. i.eq.nnt) ddc1(j)=c(j,1) if (nct.lt.nres .and. i.eq.nct) ddcn(j)=c(j,nres) ddc(j)=c(j,i) ddx(j)=c(j,i+nres) dcnorm_safe1(j)=dc_norm(j,i-1) dcnorm_safe2(j)=dc_norm(j,i) dxnorm_safe(j)=dc_norm(j,i+nres) enddo do j=1,3 c(j,i)=ddc(j)+aincr if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=c(j,1)+aincr if (nct.lt.nres .and. i.eq.nct) c(j,nres)=c(j,nres)+aincr if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1) dc(j,i)=c(j,i+1)-c(j,i) dc(j,i+nres)=c(j,i+nres)-c(j,i) call int_from_cart1(.false.) if (.not.split_ene) then call etotal(energia1) etot1=energia1(0) c write (iout,*) "ij",i,j," etot1",etot1 else !- split gradient call etotal_long(energia1) etot11=energia1(0) call etotal_short(energia1) etot12=energia1(0) endif !- end split gradient ! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1 c(j,i)=ddc(j)-aincr if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j)-aincr if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j)-aincr if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1) dc(j,i)=c(j,i+1)-c(j,i) dc(j,i+nres)=c(j,i+nres)-c(j,i) call int_from_cart1(.false.) if (.not.split_ene) then call etotal(energia1) etot2=energia1(0) c write (iout,*) "ij",i,j," etot2",etot2 ggg(j)=(etot1-etot2)/(2*aincr) else !- split gradient call etotal_long(energia1) etot21=energia1(0) ggg(j)=(etot11-etot21)/(2*aincr) call etotal_short(energia1) etot22=energia1(0) ggg1(j)=(etot12-etot22)/(2*aincr) !- end split gradient ! write (iout,*) "etot21",etot21," etot22",etot22 endif ! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2 c(j,i)=ddc(j) if (nnt.gt.1 .and. i.eq.nnt) c(j,1)=ddc1(j) if (nct.lt.nres .and. i.eq.nct) c(j,nres)=ddcn(j) if (i.gt.1) dc(j,i-1)=c(j,i)-c(j,i-1) dc(j,i)=c(j,i+1)-c(j,i) dc(j,i+nres)=c(j,i+nres)-c(j,i) dc_norm(j,i-1)=dcnorm_safe1(j) dc_norm(j,i)=dcnorm_safe2(j) dc_norm(j,i+nres)=dxnorm_safe(j) enddo do j=1,3 c(j,i+nres)=ddx(j)+aincr dc(j,i+nres)=c(j,i+nres)-c(j,i) call int_from_cart1(.false.) if (.not.split_ene) then call etotal(energia1) etot1=energia1(0) else !- split gradient call etotal_long(energia1) etot11=energia1(0) call etotal_short(energia1) etot12=energia1(0) endif !- end split gradient c(j,i+nres)=ddx(j)-aincr dc(j,i+nres)=c(j,i+nres)-c(j,i) call int_from_cart1(.false.) if (.not.split_ene) then call etotal(energia1) etot2=energia1(0) ggg(j+3)=(etot1-etot2)/(2*aincr) else !- split gradient call etotal_long(energia1) etot21=energia1(0) ggg(j+3)=(etot11-etot21)/(2*aincr) call etotal_short(energia1) etot22=energia1(0) ggg1(j+3)=(etot12-etot22)/(2*aincr) !- end split gradient endif ! write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot2",etot2 c(j,i+nres)=ddx(j) dc(j,i+nres)=c(j,i+nres)-c(j,i) dc_norm(j,i+nres)=dxnorm_safe(j) call int_from_cart1(.false.) 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 #else 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 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) write (iout,*) "enter cartgrad" call flush(iout) call cartgrad write (iout,*) "exit cartgrad" call flush(iout) icall =1 write (iout,*) "longrange grad" do i=1,nres write (iout,'(i5,3f10.5,5x,3f10.5)') 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) write (iout,*) "enter cartgrad" call flush(iout) call cartgrad write (iout,*) "exit cartgrad" call flush(iout) icall =1 write (iout,*) "shortrange grad" do i=1,nres write (iout,'(i5,3f10.5,5x,3f10.5)') 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 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 #endif 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 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 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) call gradient(nvar,x,nf,gana,uiparm,urparm,fdum) 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 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