+++ /dev/null
- subroutine check_cartgrad
-C Check the gradient of Cartesian coordinates in internal coordinates.
- implicit real*8 (a-h,o-z)
- include 'DIMENSIONS'
- include 'COMMON.IOUNITS'
- include 'COMMON.VAR'
- include 'COMMON.CHAIN'
- include 'COMMON.GEO'
- include 'COMMON.LOCAL'
- include 'COMMON.DERIV'
- dimension temp(6,maxres),xx(3),gg(3)
- indmat(i,j)=((2*(nres-2)-i)*(i-1))/2+j-1
-*
-* Check the gradient of the virtual-bond and SC vectors in the internal
-* coordinates.
-*
- aincr=1.0d-7
- aincr2=5.0d-8
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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 real*8 (a-h,o-z)
- include 'DIMENSIONS'
- include 'COMMON.CHAIN'
- include 'COMMON.DERIV'
- include 'COMMON.IOUNITS'
- include 'COMMON.VAR'
- include 'COMMON.CONTACTS'
-#ifdef MOMENT
- include 'COMMON.CONTACTS.MOMENT'
-#endif
- common /srutu/ icall
- dimension ggg(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),g(maxvar)
- dimension grad_s(6,maxres)
- double precision energia(0:n_ene),energia1(0:n_ene)
- integer uiparm(1)
- double precision urparm(1)
- external fdum
- icg=1
- nf=0
- nfl=0
- call zerograd
- aincr=1.0D-7
- print '(a)','CG processor',me,' calling CHECK_CART.'
- 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
-c----------------------------------------------------------------------------
- subroutine check_ecartint
-C Check the gradient of the energy in Cartesian coordinates.
- implicit real*8 (a-h,o-z)
- include 'DIMENSIONS'
- include 'COMMON.CONTROL'
- include 'COMMON.CHAIN'
- include 'COMMON.DERIV'
- include 'COMMON.IOUNITS'
- include 'COMMON.VAR'
- include 'COMMON.CONTACTS'
-#ifdef MOMENT
- include 'COMMON.CONTACTS.MOMENT'
-#endif
- include 'COMMON.MD_'
- include 'COMMON.LOCAL'
- include 'COMMON.SPLITELE'
- common /srutu/ icall
- dimension ggg(6),ggg1(6),cc(3),xx(3),ddc(3),ddx(3),x(maxvar),
- & g(maxvar)
- dimension dcnorm_safe(3),dxnorm_safe(3)
- dimension 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
- r_cut=2.0d0
- rlambd=0.3d0
- icg=1
- nf=0
- nfl=0
- call intout
-c call intcartderiv
-c call checkintcartgrad
- call zerograd
- aincr=1.0D-5
- write(iout,*) 'Calling CHECK_ECARTINT.'
- 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))
- call flush(iout)
- write (iout,*) "enter cartgrad"
- call flush(iout)
- call cartgrad
- 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
- 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
- write (iout,*) "split_ene not supported"
-c call zerograd
-c call etotal_long(energia(0))
-c call enerprint(energia(0))
-c call flush(iout)
-c write (iout,*) "enter cartgrad"
-c call flush(iout)
-c call cartgrad
-c write (iout,*) "exit cartgrad"
-c call flush(iout)
-c icall =1
-c write (iout,*) "longrange grad"
-c do i=1,nres
-c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
-c & (gxcart(j,i),j=1,3)
-c enddo
-c do j=1,3
-c grad_s(j,0)=gcart(j,0)
-c enddo
-c do i=1,nres
-c do j=1,3
-c grad_s(j,i)=gcart(j,i)
-c grad_s(j+3,i)=gxcart(j,i)
-c enddo
-c enddo
-c call zerograd
-c call etotal_short(energia(0))
-c call enerprint(energia(0))
-c call flush(iout)
-c write (iout,*) "enter cartgrad"
-c call flush(iout)
-c call cartgrad
-c write (iout,*) "exit cartgrad"
-c call flush(iout)
-c icall =1
-c write (iout,*) "shortrange grad"
-c do i=1,nres
-c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
-c & (gxcart(j,i),j=1,3)
-c enddo
-c do j=1,3
-c grad_s1(j,0)=gcart(j,0)
-c enddo
-c do i=1,nres
-c do j=1,3
-c grad_s1(j,i)=gcart(j,i)
-c grad_s1(j+3,i)=gxcart(j,i)
-c enddo
-c enddo
- endif
- write (iout,'(/a/)') 'Gradient in virtual-bond and SC vectors'
- do i=0,nres
- 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
-c call etotal_long(energia1(0))
-c etot11=energia1(0)
-c call etotal_short(energia1(0))
-c 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 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
-c call etotal_long(energia1(0))
-c etot21=energia1(0)
-c ggg(j)=(etot11-etot21)/(2*aincr)
-c call etotal_short(energia1(0))
-c etot22=energia1(0)
-c 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
-c call etotal_long(energia1(0))
-c etot11=energia1(0)
-c call etotal_short(energia1(0))
-c 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
-c call etotal_long(energia1(0))
-c etot21=energia1(0)
-c ggg(j+3)=(etot11-etot21)/(2*aincr)
-c call etotal_short(energia1(0))
-c etot22=energia1(0)
-c 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 real*8 (a-h,o-z)
- 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
- 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
- dnorm1=dist(i-1,i)
- dnorm2=dist(i,i+1)
- 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) phi(i+1)=beta(i-2,i-1,i,i+1)
- omeg(i)=beta(nres+i,i,maxres2,i+1)
- alph(i)=alpha(nres+i,i,maxres2)
- theta(i+1)=alpha(i-1,i,i+1)
- vbld(i)=dist(i-1,i)
- vbld_inv(i)=1.0d0/vbld(i)
- vbld(nres+i)=dist(nres+i,i)
- 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 real*8 (a-h,o-z)
- include 'DIMENSIONS'
- include 'COMMON.CHAIN'
- include 'COMMON.DERIV'
- include 'COMMON.IOUNITS'
- include 'COMMON.VAR'
- include 'COMMON.GEO'
- common /srutu/ icall
- dimension 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
- external fdum
- call zerograd
- aincr=1.0D-7
- print '(a)','Calling CHECK_INT.'
- 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
- call etotal(energia1(0))
- etot1=energia1(0)
- x(i)=xi+0.5D0*aincr
- call var_to_geom(nvar,x)
- call chainbuild
- 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
projects / unres.git / blobdiff