--- /dev/null
+ module geometry
+!-----------------------------------------------------------------------------
+ use io_units
+ use names
+ use math
+ use MPI_data
+ use geometry_data
+ use control_data
+ use energy_data
+ implicit none
+!-----------------------------------------------------------------------------
+! commom.bounds
+! common /bounds/
+!-----------------------------------------------------------------------------
+! commom.chain
+! common /chain/
+! common /rotmat/
+ real(kind=8),dimension(:,:,:),allocatable :: t,r !(3,3,maxres)
+!-----------------------------------------------------------------------------
+! common.geo
+! common /geo/
+!-----------------------------------------------------------------------------
+! common.locmove
+! Variables (set in init routine) never modified by local_move
+! common /loc_const/
+ integer :: init_called
+ logical :: locmove_output
+ real(kind=8) :: min_theta, max_theta
+ real(kind=8) :: dmin2,dmax2
+ real(kind=8) :: flag,small,small2
+! Workspace for local_move
+! common /loc_work/
+ integer :: a_n,b_n,res_n
+ real(kind=8),dimension(0:7) :: a_ang
+ real(kind=8),dimension(0:3) :: b_ang
+ real(kind=8),dimension(0:11) :: res_ang
+ logical,dimension(0:2,0:7) :: a_tab
+ logical,dimension(0:2,0:3) :: b_tab
+ logical,dimension(0:2,0:2,0:11) :: res_tab
+!-----------------------------------------------------------------------------
+! integer,dimension(:),allocatable :: itype_pdb !(maxres) initialize in molread
+!-----------------------------------------------------------------------------
+!
+!
+!-----------------------------------------------------------------------------
+ contains
+!-----------------------------------------------------------------------------
+! arcos.f
+!-----------------------------------------------------------------------------
+ real(kind=8) function ARCOS(X)
+! implicit real*8 (a-h,o-z)
+! include 'COMMON.GEO'
+!el local variables
+ real(kind=8) :: x
+ IF (DABS(X).LT.1.0D0) GOTO 1
+ ARCOS=PIPOL*(1.0d0-DSIGN(1.0D0,X))
+ RETURN
+ 1 ARCOS=DACOS(X)
+ return
+ end function ARCOS
+!-----------------------------------------------------------------------------
+! chainbuild.F
+!-----------------------------------------------------------------------------
+ subroutine chainbuild
+!
+! Build the virtual polypeptide chain. Side-chain centroids are moveable.
+! As of 2/17/95.
+!
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.LOCAL'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.NAMES'
+! include 'COMMON.INTERACT'
+ logical :: lprn
+!el local variables
+ integer :: i,j
+ real(kind=8) :: be,be1,alfai
+ integer :: nres2
+ nres2=2*nres
+! Set lprn=.true. for debugging
+ lprn = .false.
+!
+! Define the origin and orientation of the coordinate system and locate the
+! first three CA's and SC(2).
+!
+!elwrite(iout,*)"in chainbuild"
+ call orig_frame
+!elwrite(iout,*)"after orig_frame"
+!
+! Build the alpha-carbon chain.
+!
+ do i=4,nres
+ call locate_next_res(i)
+ enddo
+!elwrite(iout,*)"after locate_next_res"
+!
+! First and last SC must coincide with the corresponding CA.
+!
+ do j=1,3
+ dc(j,nres+1)=0.0D0
+ dc_norm(j,nres+1)=0.0D0
+ dc(j,nres+nres)=0.0D0
+ dc_norm(j,nres+nres)=0.0D0
+ c(j,nres+1)=c(j,1)
+ c(j,nres+nres)=c(j,nres)
+ enddo
+!
+! Temporary diagnosis
+!
+ if (lprn) then
+
+ call cartprint
+ write (iout,'(/a)') 'Recalculated internal coordinates'
+ do i=2,nres-1
+ do j=1,3
+ c(j,nres2+2)=0.5D0*(c(j,i-1)+c(j,i+1)) !maxres2=2*maxres
+ enddo
+ be=0.0D0
+ if (i.gt.3) be=rad2deg*beta(i-3,i-2,i-1,i)
+ be1=rad2deg*beta(nres+i,i,nres2+2,i+1)
+ alfai=0.0D0
+ if (i.gt.2) alfai=rad2deg*alpha(i-2,i-1,i)
+ write (iout,1212) restyp(itype(i)),i,dist(i-1,i),&
+ alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,nres2+2),be1
+ enddo
+ 1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
+
+ endif
+
+ return
+ end subroutine chainbuild
+!-----------------------------------------------------------------------------
+ subroutine orig_frame
+!
+! Define the origin and orientation of the coordinate system and locate
+! the first three atoms.
+!
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.LOCAL'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+!el local variables
+ integer :: i,j
+ real(kind=8) :: cost,sint
+
+!el allocate(t(3,3,nres)) !(3,3,maxres)
+!el allocate(r(3,3,nres)) !(3,3,maxres)
+!el allocate(rt(3,3,nres)) !(3,3,maxres)
+!el allocate(dc_norm(3,0:2*nres)) !(3,0:maxres2)
+!el allocate(prod(3,3,nres)) !(3,3,maxres)
+
+ cost=dcos(theta(3))
+ sint=dsin(theta(3))
+ t(1,1,1)=-cost
+ t(1,2,1)=-sint
+ t(1,3,1)= 0.0D0
+ t(2,1,1)=-sint
+ t(2,2,1)= cost
+ t(2,3,1)= 0.0D0
+ t(3,1,1)= 0.0D0
+ t(3,2,1)= 0.0D0
+ t(3,3,1)= 1.0D0
+ r(1,1,1)= 1.0D0
+ r(1,2,1)= 0.0D0
+ r(1,3,1)= 0.0D0
+ r(2,1,1)= 0.0D0
+ r(2,2,1)= 1.0D0
+ r(2,3,1)= 0.0D0
+ r(3,1,1)= 0.0D0
+ r(3,2,1)= 0.0D0
+ r(3,3,1)= 1.0D0
+ do i=1,3
+ do j=1,3
+ rt(i,j,1)=t(i,j,1)
+ enddo
+ enddo
+ do i=1,3
+ do j=1,3
+ prod(i,j,1)=0.0D0
+ prod(i,j,2)=t(i,j,1)
+ enddo
+ prod(i,i,1)=1.0D0
+ enddo
+ c(1,1)=0.0D0
+ c(2,1)=0.0D0
+ c(3,1)=0.0D0
+ c(1,2)=vbld(2)
+ c(2,2)=0.0D0
+ c(3,2)=0.0D0
+ dc(1,0)=0.0d0
+ dc(2,0)=0.0D0
+ dc(3,0)=0.0D0
+ dc(1,1)=vbld(2)
+ dc(2,1)=0.0D0
+ dc(3,1)=0.0D0
+ dc_norm(1,0)=0.0D0
+ dc_norm(2,0)=0.0D0
+ dc_norm(3,0)=0.0D0
+ dc_norm(1,1)=1.0D0
+ dc_norm(2,1)=0.0D0
+ dc_norm(3,1)=0.0D0
+ do j=1,3
+ dc_norm(j,2)=prod(j,1,2)
+ dc(j,2)=vbld(3)*prod(j,1,2)
+ c(j,3)=c(j,2)+dc(j,2)
+ enddo
+ call locate_side_chain(2)
+ return
+ end subroutine orig_frame
+!-----------------------------------------------------------------------------
+ subroutine locate_next_res(i)
+!
+! Locate CA(i) and SC(i-1)
+!
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.LOCAL'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.NAMES'
+! include 'COMMON.INTERACT'
+!
+! Define the rotation matrices corresponding to CA(i)
+!
+!el local variables
+ integer :: i,j
+ real(kind=8) :: theti,phii
+ real(kind=8) :: cost,sint,cosphi,sinphi
+#ifdef OSF
+#ifdef WHAM_RUN
+ theti=theta(i)
+ icrc=0
+ call proc_proc(theti,icrc)
+ if(icrc.eq.1)theti=100.0
+ phii=phi(i)
+ icrc=0
+ call proc_proc(phii,icrc)
+ if(icrc.eq.1)phii=180.0
+#else
+ theti=theta(i)
+ if (theti.ne.theti) theti=100.0
+ phii=phi(i)
+ if (phii.ne.phii) phii=180.0
+#endif
+#else
+ theti=theta(i)
+ phii=phi(i)
+#endif
+ cost=dcos(theti)
+ sint=dsin(theti)
+ cosphi=dcos(phii)
+ sinphi=dsin(phii)
+! Define the matrices of the rotation about the virtual-bond valence angles
+! theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
+! program), R(i,j,k), and, the cumulative matrices of rotation RT
+ t(1,1,i-2)=-cost
+ t(1,2,i-2)=-sint
+ t(1,3,i-2)= 0.0D0
+ t(2,1,i-2)=-sint
+ t(2,2,i-2)= cost
+ t(2,3,i-2)= 0.0D0
+ t(3,1,i-2)= 0.0D0
+ t(3,2,i-2)= 0.0D0
+ t(3,3,i-2)= 1.0D0
+ r(1,1,i-2)= 1.0D0
+ r(1,2,i-2)= 0.0D0
+ r(1,3,i-2)= 0.0D0
+ r(2,1,i-2)= 0.0D0
+ r(2,2,i-2)=-cosphi
+ r(2,3,i-2)= sinphi
+ r(3,1,i-2)= 0.0D0
+ r(3,2,i-2)= sinphi
+ r(3,3,i-2)= cosphi
+ rt(1,1,i-2)=-cost
+ rt(1,2,i-2)=-sint
+ rt(1,3,i-2)=0.0D0
+ rt(2,1,i-2)=sint*cosphi
+ rt(2,2,i-2)=-cost*cosphi
+ rt(2,3,i-2)=sinphi
+ rt(3,1,i-2)=-sint*sinphi
+ rt(3,2,i-2)=cost*sinphi
+ rt(3,3,i-2)=cosphi
+ call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
+ do j=1,3
+ dc_norm(j,i-1)=prod(j,1,i-1)
+ dc(j,i-1)=vbld(i)*prod(j,1,i-1)
+ c(j,i)=c(j,i-1)+dc(j,i-1)
+ enddo
+!d print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3)
+!
+! Now calculate the coordinates of SC(i-1)
+!
+ call locate_side_chain(i-1)
+ return
+ end subroutine locate_next_res
+!-----------------------------------------------------------------------------
+ subroutine locate_side_chain(i)
+!
+! Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
+!
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.LOCAL'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.NAMES'
+! include 'COMMON.INTERACT'
+ integer :: i,j,k
+ real(kind=8),dimension(3) :: xx
+ real(kind=8) :: alphi,omegi,theta2
+ real(kind=8) :: dsci,dsci_inv,sinalphi,cosalphi,cosomegi,sinomegi
+ real(kind=8) :: xp,yp,zp,cost2,sint2,rj
+! dsci=dsc(itype(i))
+! dsci_inv=dsc_inv(itype(i))
+ dsci=vbld(i+nres)
+ dsci_inv=vbld_inv(i+nres)
+#ifdef OSF
+ alphi=alph(i)
+ omegi=omeg(i)
+#ifdef WHAM_RUN
+! detecting NaNQ
+ icrc=0
+ call proc_proc(alphi,icrc)
+ if(icrc.eq.1)alphi=100.0
+ icrc=0
+ call proc_proc(omegi,icrc)
+ if(icrc.eq.1)omegi=-100.0
+#else
+ if (alphi.ne.alphi) alphi=100.0
+ if (omegi.ne.omegi) omegi=-100.0
+#endif
+#else
+ alphi=alph(i)
+ omegi=omeg(i)
+#endif
+ cosalphi=dcos(alphi)
+ sinalphi=dsin(alphi)
+ cosomegi=dcos(omegi)
+ sinomegi=dsin(omegi)
+ xp= dsci*cosalphi
+ yp= dsci*sinalphi*cosomegi
+ zp=-dsci*sinalphi*sinomegi
+! Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its
+! X-axis aligned with the vector DC(*,i)
+ theta2=pi-0.5D0*theta(i+1)
+ cost2=dcos(theta2)
+ sint2=dsin(theta2)
+ xx(1)= xp*cost2+yp*sint2
+ xx(2)=-xp*sint2+yp*cost2
+ xx(3)= zp
+!d print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
+!d & xp,yp,zp,(xx(k),k=1,3)
+ do j=1,3
+ xloc(j,i)=xx(j)
+ enddo
+! Bring the SC vectors to the common coordinate system.
+ xx(1)=xloc(1,i)
+ xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1)
+ xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1)
+ do j=1,3
+ xrot(j,i)=xx(j)
+ enddo
+ do j=1,3
+ rj=0.0D0
+ do k=1,3
+ rj=rj+prod(j,k,i-1)*xx(k)
+ enddo
+ dc(j,nres+i)=rj
+ dc_norm(j,nres+i)=rj*dsci_inv
+ c(j,nres+i)=c(j,i)+rj
+ enddo
+ return
+ end subroutine locate_side_chain
+!-----------------------------------------------------------------------------
+! checkder_p.F
+!-----------------------------------------------------------------------------
+ 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
+!el local variables
+ integer :: i,j
+ real(kind=8) :: dnorm1,dnorm2,be
+ integer :: nres2
+ nres2=2*nres
+ if (lprn) write (iout,'(/a)') 'Recalculated internal coordinates'
+#ifdef TIMING
+ time01=MPI_Wtime()
+#endif
+
+#ifdef WHAM_RUN
+ vbld(nres+1)=0.0d0
+!write(iout,*)"geometry warring, vbld=",(vbld(i),i=1,nres+1)
+ vbld(2*nres)=0.0d0
+ vbld_inv(nres+1)=0.0d0
+ vbld_inv(2*nres)=0.0d0
+#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,nres2+2)=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,nres2+2,i+1)
+ alph(i)=alpha(nres+i,i,nres2+2)
+ 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
+!d write(iout,*) "iint_start",iint_start," iint_count",
+!d & (iint_count(i),i=0,nfgtasks-1)," iint_displ",
+!d & (iint_displ(i),i=0,nfgtasks-1)
+!d write (iout,*) "Gather vbld backbone"
+!d 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)
+!d write (iout,*) "Gather vbld_inv"
+!d 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)
+!d write (iout,*) "Gather vbld side chain"
+!d 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)
+!d write (iout,*) "Gather vbld_inv side chain"
+!d 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)
+!d write (iout,*) "Gather theta"
+!d 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)
+!d write (iout,*) "Gather phi"
+!d 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
+!d write (iout,*) "Gather alph"
+!d 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)
+!d write (iout,*) "Gather omeg"
+!d 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
+!#ifdef WHAM_RUN
+#if defined(WHAM_RUN) || defined(CLUSTER)
+ dc(j,i)=c(j,i+1)-c(j,i)
+#endif
+ dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+ enddo
+ enddo
+ do i=2,nres-1
+ do j=1,3
+!#ifdef WHAM_RUN
+#if defined(WHAM_RUN) || defined(CLUSTER)
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+#endif
+ 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 subroutine int_from_cart1
+#if !defined(WHAM_RUN) && !defined(CLUSTER)
+!-----------------------------------------------------------------------------
+! check_sc_distr.f
+!-----------------------------------------------------------------------------
+ subroutine check_sc_distr
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.TIME1'
+! include 'COMMON.INTERACT'
+! include 'COMMON.NAMES'
+! include 'COMMON.GEO'
+! include 'COMMON.HEADER'
+! include 'COMMON.CONTROL'
+ logical :: fail
+ real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
+ real(kind=8) :: hrtime,mintime,sectime
+ integer,parameter :: MaxSample=10000000
+ real(kind=8),parameter :: delt=1.0D0/MaxSample
+ real(kind=8),dimension(0:72,0:90) :: prob
+!el local variables
+ integer :: it,i,j,isample,indal,indom
+ real(kind=8) :: al,om,dV
+ dV=2.0D0*5.0D0*deg2rad*deg2rad
+ print *,'dv=',dv
+ do 10 it=1,1
+ if (it.eq.10) goto 10
+ open (20,file=restyp(it)//'_distr.sdc',status='unknown')
+ call gen_side(it,90.0D0 * deg2rad,al,om,fail)
+ close (20)
+ goto 10
+ open (20,file=restyp(it)//'_distr1.sdc',status='unknown')
+ do i=0,90
+ do j=0,72
+ prob(j,i)=0.0D0
+ enddo
+ enddo
+ do isample=1,MaxSample
+ call gen_side(it,90.0D0 * deg2rad,al,om,fail)
+ indal=rad2deg*al/2
+ indom=(rad2deg*om+180.0D0)/5
+ prob(indom,indal)=prob(indom,indal)+delt
+ enddo
+ do i=45,90
+ do j=0,72
+ write (20,'(2f10.3,1pd15.5)') 2*i+0.0D0,5*j-180.0D0,&
+ prob(j,i)/dV
+ enddo
+ enddo
+ 10 continue
+ return
+ end subroutine check_sc_distr
+#endif
+!-----------------------------------------------------------------------------
+! convert.f
+!-----------------------------------------------------------------------------
+ subroutine geom_to_var(n,x)
+!
+! Transfer the geometry parameters to the variable array.
+! The positions of variables are as follows:
+! 1. Virtual-bond torsional angles: 1 thru nres-3
+! 2. Virtual-bond valence angles: nres-2 thru 2*nres-5
+! 3. The polar angles alpha of local SC orientation: 2*nres-4 thru
+! 2*nres-4+nside
+! 4. The torsional angles omega of SC orientation: 2*nres-4+nside+1
+! thru 2*nre-4+2*nside
+!
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.VAR'
+! include 'COMMON.GEO'
+! include 'COMMON.CHAIN'
+ integer :: n,i
+ real(kind=8),dimension(n) :: x
+!d print *,'nres',nres,' nphi',nphi,' ntheta',ntheta,' nvar',nvar
+ do i=4,nres
+ x(i-3)=phi(i)
+!d print *,i,i-3,phi(i)
+ enddo
+ if (n.eq.nphi) return
+ do i=3,nres
+ x(i-2+nphi)=theta(i)
+!d print *,i,i-2+nphi,theta(i)
+ enddo
+ if (n.eq.nphi+ntheta) return
+ do i=2,nres-1
+ if (ialph(i,1).gt.0) then
+ x(ialph(i,1))=alph(i)
+ x(ialph(i,1)+nside)=omeg(i)
+!d print *,i,ialph(i,1),ialph(i,1)+nside,alph(i),omeg(i)
+ endif
+ enddo
+ return
+ end subroutine geom_to_var
+!-----------------------------------------------------------------------------
+ subroutine var_to_geom(n,x)
+!
+! Update geometry parameters according to the variable array.
+!
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.VAR'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! include 'COMMON.IOUNITS'
+ integer :: n,i,ii
+ real(kind=8),dimension(n) :: x
+ logical :: change !,reduce
+!el alph=0.0d0
+!el omeg=0.0d0
+!el phi=0.0d0
+!el theta=0.0d0
+
+ change=reduce(x)
+ if (n.gt.nphi+ntheta) then
+ do i=1,nside
+ ii=ialph(i,2)
+ alph(ii)=x(nphi+ntheta+i)
+ omeg(ii)=pinorm(x(nphi+ntheta+nside+i))
+!elwrite(iout,*) "alph",ii,alph
+!elwrite(iout,*) "omeg",ii,omeg
+ enddo
+ endif
+ do i=4,nres
+ phi(i)=x(i-3)
+!elwrite(iout,*) "phi",i,phi
+ enddo
+ if (n.eq.nphi) return
+ do i=3,nres
+ theta(i)=x(i-2+nphi)
+!elwrite(iout,*) "theta",i,theta
+ if (theta(i).eq.pi) theta(i)=0.99d0*pi
+ x(i-2+nphi)=theta(i)
+ enddo
+ return
+ end subroutine var_to_geom
+!-----------------------------------------------------------------------------
+ logical function convert_side(alphi,omegi)
+! implicit none
+ real(kind=8) :: alphi,omegi
+!el real(kind=8) :: pinorm
+! include 'COMMON.GEO'
+ convert_side=.false.
+! Apply periodicity restrictions.
+ if (alphi.gt.pi) then
+ alphi=dwapi-alphi
+ omegi=pinorm(omegi+pi)
+ convert_side=.true.
+ endif
+ return
+ end function convert_side
+!-----------------------------------------------------------------------------
+ logical function reduce(x)
+!
+! Apply periodic restrictions to variables.
+!
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.VAR'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+ logical :: zm,zmiana !,convert_side
+ real(kind=8),dimension(nvar) :: x
+ integer :: i,ii,iii
+ zmiana=.false.
+ do i=4,nres
+ x(i-3)=pinorm(x(i-3))
+ enddo
+ if (nvar.gt.nphi+ntheta) then
+ do i=1,nside
+ ii=nphi+ntheta+i
+ iii=ii+nside
+ x(ii)=thetnorm(x(ii))
+ x(iii)=pinorm(x(iii))
+! Apply periodic restrictions.
+ zm=convert_side(x(ii),x(iii))
+ zmiana=zmiana.or.zm
+ enddo
+ endif
+ if (nvar.eq.nphi) return
+ do i=3,nres
+ ii=i-2+nphi
+ iii=i-3
+ x(ii)=dmod(x(ii),dwapi)
+! Apply periodic restrictions.
+ if (x(ii).gt.pi) then
+ zmiana=.true.
+ x(ii)=dwapi-x(ii)
+ if (iii.gt.0) x(iii)=pinorm(x(iii)+pi)
+ if (i.lt.nres) x(iii+1)=pinorm(x(iii+1)+pi)
+ ii=ialph(i-1,1)
+ if (ii.gt.0) then
+ x(ii)=dmod(pi-x(ii),dwapi)
+ x(ii+nside)=pinorm(-x(ii+nside))
+ zm=convert_side(x(ii),x(ii+nside))
+ endif
+ else if (x(ii).lt.-pi) then
+ zmiana=.true.
+ x(ii)=dwapi+x(ii)
+ ii=ialph(i-1,1)
+ if (ii.gt.0) then
+ x(ii)=dmod(pi-x(ii),dwapi)
+ x(ii+nside)=pinorm(-pi-x(ii+nside))
+ zm=convert_side(x(ii),x(ii+nside))
+ endif
+ else if (x(ii).lt.0.0d0) then
+ zmiana=.true.
+ x(ii)=-x(ii)
+ if (iii.gt.0) x(iii)=pinorm(x(iii)+pi)
+ if (i.lt.nres) x(iii+1)=pinorm(x(iii+1)+pi)
+ ii=ialph(i-1,1)
+ if (ii.gt.0) then
+ x(ii+nside)=pinorm(-x(ii+nside))
+ zm=convert_side(x(ii),x(ii+nside))
+ endif
+ endif
+ enddo
+ reduce=zmiana
+ return
+ end function reduce
+!-----------------------------------------------------------------------------
+ real(kind=8) function thetnorm(x)
+! This function puts x within [0,2Pi].
+ implicit none
+ real(kind=8) :: x,xx
+! include 'COMMON.GEO'
+ xx=dmod(x,dwapi)
+ if (xx.lt.0.0d0) xx=xx+dwapi
+ if (xx.gt.0.9999d0*pi) xx=0.9999d0*pi
+ thetnorm=xx
+ return
+ end function thetnorm
+#if !defined(WHAM_RUN) && !defined(CLUSTER)
+!-----------------------------------------------------------------------------
+ subroutine var_to_geom_restr(n,xx)
+!
+! Update geometry parameters according to the variable array.
+!
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.VAR'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! include 'COMMON.IOUNITS'
+ integer :: n,i,ii
+ real(kind=8),dimension(6*nres) :: x,xx !(maxvar) (maxvar=6*maxres)
+ logical :: change !,reduce
+
+ call xx2x(x,xx)
+ change=reduce(x)
+ do i=1,nside
+ ii=ialph(i,2)
+ alph(ii)=x(nphi+ntheta+i)
+ omeg(ii)=pinorm(x(nphi+ntheta+nside+i))
+ enddo
+ do i=4,nres
+ phi(i)=x(i-3)
+ enddo
+ do i=3,nres
+ theta(i)=x(i-2+nphi)
+ if (theta(i).eq.pi) theta(i)=0.99d0*pi
+ x(i-2+nphi)=theta(i)
+ enddo
+ return
+ end subroutine var_to_geom_restr
+!-----------------------------------------------------------------------------
+! gen_rand_conf.F
+!-----------------------------------------------------------------------------
+ subroutine gen_rand_conf(nstart,*)
+! Generate random conformation or chain cut and regrowth.
+ use mcm_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.LOCAL'
+! include 'COMMON.VAR'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.MCM'
+! include 'COMMON.GEO'
+! include 'COMMON.CONTROL'
+ logical :: back,fail !overlap,
+!el local variables
+ integer :: i,nstart,maxsi,nsi,maxnit,nit,niter
+ integer :: it1,it2,it,j
+!d print *,' CG Processor',me,' maxgen=',maxgen
+ maxsi=100
+!d write (iout,*) 'Gen_Rand_conf: nstart=',nstart
+ if (nstart.lt.5) then
+ it1=iabs(itype(2))
+ phi(4)=gen_phi(4,iabs(itype(2)),iabs(itype(3)))
+! write(iout,*)'phi(4)=',rad2deg*phi(4)
+ if (nstart.lt.3) theta(3)=gen_theta(iabs(itype(2)),pi,phi(4))
+! write(iout,*)'theta(3)=',rad2deg*theta(3)
+ if (it1.ne.10) then
+ nsi=0
+ fail=.true.
+ do while (fail.and.nsi.le.maxsi)
+ call gen_side(it1,theta(3),alph(2),omeg(2),fail)
+ nsi=nsi+1
+ enddo
+ if (nsi.gt.maxsi) return 1
+ endif ! it1.ne.10
+ call orig_frame
+ i=4
+ nstart=4
+ else
+ i=nstart
+ nstart=max0(i,4)
+ endif
+
+ maxnit=0
+
+ nit=0
+ niter=0
+ back=.false.
+ do while (i.le.nres .and. niter.lt.maxgen)
+ if (i.lt.nstart) then
+ if(iprint.gt.1) then
+ write (iout,'(/80(1h*)/2a/80(1h*))') &
+ 'Generation procedure went down to ',&
+ 'chain beginning. Cannot continue...'
+ write (*,'(/80(1h*)/2a/80(1h*))') &
+ 'Generation procedure went down to ',&
+ 'chain beginning. Cannot continue...'
+ endif
+ return 1
+ endif
+ it1=iabs(itype(i-1))
+ it2=iabs(itype(i-2))
+ it=iabs(itype(i))
+! print *,'Gen_Rand_Conf: i=',i,' it=',it,' it1=',it1,' it2=',it2,
+! & ' nit=',nit,' niter=',niter,' maxgen=',maxgen
+ phi(i+1)=gen_phi(i+1,it1,it)
+ if (back) then
+ phi(i)=gen_phi(i+1,it2,it1)
+! print *,'phi(',i,')=',phi(i)
+ theta(i-1)=gen_theta(it2,phi(i-1),phi(i))
+ if (it2.ne.10) then
+ nsi=0
+ fail=.true.
+ do while (fail.and.nsi.le.maxsi)
+ call gen_side(it2,theta(i-1),alph(i-2),omeg(i-2),fail)
+ nsi=nsi+1
+ enddo
+ if (nsi.gt.maxsi) return 1
+ endif
+ call locate_next_res(i-1)
+ endif
+ theta(i)=gen_theta(it1,phi(i),phi(i+1))
+ if (it1.ne.10) then
+ nsi=0
+ fail=.true.
+ do while (fail.and.nsi.le.maxsi)
+ call gen_side(it1,theta(i),alph(i-1),omeg(i-1),fail)
+ nsi=nsi+1
+ enddo
+ if (nsi.gt.maxsi) return 1
+ endif
+ call locate_next_res(i)
+ if (overlap(i-1)) then
+ if (nit.lt.maxnit) then
+ back=.true.
+ nit=nit+1
+ else
+ nit=0
+ if (i.gt.3) then
+ back=.true.
+ i=i-1
+ else
+ write (iout,'(a)') &
+ 'Cannot generate non-overlaping conformation. Increase MAXNIT.'
+ write (*,'(a)') &
+ 'Cannot generate non-overlaping conformation. Increase MAXNIT.'
+ return 1
+ endif
+ endif
+ else
+ back=.false.
+ nit=0
+ i=i+1
+ endif
+ niter=niter+1
+ enddo
+ if (niter.ge.maxgen) then
+ write (iout,'(a,2i5)') &
+ 'Too many trials in conformation generation',niter,maxgen
+ write (*,'(a,2i5)') &
+ 'Too many trials in conformation generation',niter,maxgen
+ return 1
+ endif
+ do j=1,3
+ c(j,nres+1)=c(j,1)
+ c(j,nres+nres)=c(j,nres)
+ enddo
+ return
+ end subroutine gen_rand_conf
+!-----------------------------------------------------------------------------
+ logical function overlap(i)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.FFIELD'
+ integer :: i,j,iti,itj,iteli,itelj,k
+ real(kind=8) :: redfac,rcomp
+ integer :: nres2
+ nres2=2*nres
+ data redfac /0.5D0/
+ overlap=.false.
+ iti=iabs(itype(i))
+ if (iti.gt.ntyp) return
+! Check for SC-SC overlaps.
+!d print *,'nnt=',nnt,' nct=',nct
+ do j=nnt,i-1
+ itj=iabs(itype(j))
+ if (j.lt.i-1 .or. ipot.ne.4) then
+ rcomp=sigmaii(iti,itj)
+ else
+ rcomp=sigma(iti,itj)
+ endif
+!d print *,'j=',j
+ if (dist(nres+i,nres+j).lt.redfac*rcomp) then
+ overlap=.true.
+! print *,'overlap, SC-SC: i=',i,' j=',j,
+! & ' dist=',dist(nres+i,nres+j),' rcomp=',
+! & rcomp
+ return
+ endif
+ enddo
+! Check for overlaps between the added peptide group and the preceding
+! SCs.
+ iteli=itel(i)
+ do j=1,3
+! c(j,nres2+1)=0.5D0*(c(j,i)+c(j,i+1))
+ c(j,nres2+3)=0.5D0*(c(j,i)+c(j,i+1))
+ enddo
+ do j=nnt,i-2
+ itj=iabs(itype(j))
+!d print *,'overlap, p-Sc: i=',i,' j=',j,
+!d & ' dist=',dist(nres+j,maxres2+1)
+ if (dist(nres+j,nres2+3).lt.4.0D0*redfac) then
+ overlap=.true.
+ return
+ endif
+ enddo
+! Check for overlaps between the added side chain and the preceding peptide
+! groups.
+ do j=1,nnt-2
+ do k=1,3
+ c(k,nres2+3)=0.5D0*(c(k,j)+c(k,j+1))
+ enddo
+!d print *,'overlap, SC-p: i=',i,' j=',j,
+!d & ' dist=',dist(nres+i,maxres2+1)
+ if (dist(nres+i,nres2+3).lt.4.0D0*redfac) then
+ overlap=.true.
+ return
+ endif
+ enddo
+! Check for p-p overlaps
+ do j=1,3
+ c(j,nres2+4)=0.5D0*(c(j,i)+c(j,i+1))
+ enddo
+ do j=nnt,i-2
+ itelj=itel(j)
+ do k=1,3
+ c(k,nres2+4)=0.5D0*(c(k,j)+c(k,j+1))
+ enddo
+!d print *,'overlap, p-p: i=',i,' j=',j,
+!d & ' dist=',dist(maxres2+1,maxres2+2)
+ if(iteli.ne.0.and.itelj.ne.0)then
+ if (dist(nres2+3,nres2+4).lt.rpp(iteli,itelj)*redfac) then
+ overlap=.true.
+ return
+ endif
+ endif
+ enddo
+ return
+ end function overlap
+!-----------------------------------------------------------------------------
+ real(kind=8) function gen_phi(i,it1,it2)
+ use random, only:ran_number
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.BOUNDS'
+ integer :: i,it1,it2
+! gen_phi=ran_number(-pi,pi)
+! 8/13/98 Generate phi using pre-defined boundaries
+ gen_phi=ran_number(phibound(1,i),phibound(2,i))
+ return
+ end function gen_phi
+!-----------------------------------------------------------------------------
+ real(kind=8) function gen_theta(it,gama,gama1)
+ use random,only:binorm
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.LOCAL'
+! include 'COMMON.GEO'
+ real(kind=8),dimension(2) :: y,z
+ real(kind=8) :: theta_max,theta_min,sig,ak
+!el local variables
+ integer :: j,it,k
+ real(kind=8) :: gama,gama1,thet_pred_mean,theta_temp
+! print *,'gen_theta: it=',it
+ theta_min=0.05D0*pi
+ theta_max=0.95D0*pi
+ if (dabs(gama).gt.dwapi) then
+ y(1)=dcos(gama)
+ y(2)=dsin(gama)
+ else
+ y(1)=0.0D0
+ y(2)=0.0D0
+ endif
+ if (dabs(gama1).gt.dwapi) then
+ z(1)=dcos(gama1)
+ z(2)=dsin(gama1)
+ else
+ z(1)=0.0D0
+ z(2)=0.0D0
+ endif
+ thet_pred_mean=a0thet(it)
+ do k=1,2
+ thet_pred_mean=thet_pred_mean+athet(k,it,1,1)*y(k) &
+ +bthet(k,it,1,1)*z(k)
+ enddo
+ sig=polthet(3,it)
+ do j=2,0,-1
+ sig=sig*thet_pred_mean+polthet(j,it)
+ enddo
+ sig=0.5D0/(sig*sig+sigc0(it))
+ ak=dexp(gthet(1,it)- &
+ 0.5D0*((gthet(2,it)-thet_pred_mean)/gthet(3,it))**2)
+! print '(i5,5(1pe14.4))',it,(gthet(j,it),j=1,3)
+! print '(5(1pe14.4))',thet_pred_mean,theta0(it),sig,sig0(it),ak
+ theta_temp=binorm(thet_pred_mean,theta0(it),sig,sig0(it),ak)
+ if (theta_temp.lt.theta_min) theta_temp=theta_min
+ if (theta_temp.gt.theta_max) theta_temp=theta_max
+ gen_theta=theta_temp
+! print '(a)','Exiting GENTHETA.'
+ return
+ end function gen_theta
+!-----------------------------------------------------------------------------
+ subroutine gen_side(it,the,al,om,fail)
+ use random, only:ran_number,mult_norm1
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.LOCAL'
+! include 'COMMON.SETUP'
+! include 'COMMON.IOUNITS'
+ real(kind=8) :: MaxBoxLen=10.0D0
+ real(kind=8),dimension(3,3) :: Ap_inv,a,vec
+ real(kind=8),dimension(:,:),allocatable :: z !(3,maxlob)
+ real(kind=8),dimension(:),allocatable :: W1,detAp !(maxlob)
+ real(kind=8),dimension(:),allocatable :: sumW !(0:maxlob)
+ real(kind=8),dimension(2) :: y,cm,eig
+ real(kind=8),dimension(2,2) :: box
+ real(kind=8),dimension(100) :: work
+ real(kind=8) :: eig_limit=1.0D-8
+ real(kind=8) :: Big=10.0D0
+ logical :: lprint,fail,lcheck
+!el local variables
+ integer :: it,i,j,k,l,nlobit,ial,iom,iii,ilob
+ real(kind=8) :: the,al,om,detApi,wart,y2,wykl,radmax
+ real(kind=8) :: tant,zz1,W1i,radius,zk,fac,dV,sum,sum1
+ real(kind=8) :: which_lobe
+ lcheck=.false.
+ lprint=.false.
+ fail=.false.
+ if (the.eq.0.0D0 .or. the.eq.pi) then
+#ifdef MPI
+ write (*,'(a,i4,a,i3,a,1pe14.5)') &
+ 'CG Processor:',me,' Error in GenSide: it=',it,' theta=',the
+#else
+!d write (iout,'(a,i3,a,1pe14.5)')
+!d & 'Error in GenSide: it=',it,' theta=',the
+#endif
+ fail=.true.
+ return
+ endif
+ tant=dtan(the-pipol)
+ nlobit=nlob(it)
+ allocate(z(3,nlobit))
+ allocate(W1(nlobit))
+ allocate(detAp(nlobit))
+ allocate(sumW(0:nlobit))
+ if (lprint) then
+#ifdef MPI
+ print '(a,i4,a)','CG Processor:',me,' Enter Gen_Side.'
+ write (iout,'(a,i4,a)') 'Processor:',me,' Enter Gen_Side.'
+#endif
+ print *,'it=',it,' nlobit=',nlobit,' the=',the,' tant=',tant
+ write (iout,*) 'it=',it,' nlobit=',nlobit,' the=',the,&
+ ' tant=',tant
+ endif
+ do i=1,nlobit
+ zz1=tant-censc(1,i,it)
+ do k=1,3
+ do l=1,3
+ a(k,l)=gaussc(k,l,i,it)
+ enddo
+ enddo
+ detApi=a(2,2)*a(3,3)-a(2,3)**2
+ Ap_inv(2,2)=a(3,3)/detApi
+ Ap_inv(2,3)=-a(2,3)/detApi
+ Ap_inv(3,2)=Ap_inv(2,3)
+ Ap_inv(3,3)=a(2,2)/detApi
+ if (lprint) then
+ write (*,'(/a,i2/)') 'Cluster #',i
+ write (*,'(3(1pe14.5),5x,1pe14.5)') &
+ ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
+ write (iout,'(/a,i2/)') 'Cluster #',i
+ write (iout,'(3(1pe14.5),5x,1pe14.5)') &
+ ((a(l,k),l=1,3),censc(k,i,it),k=1,3)
+ endif
+ W1i=0.0D0
+ do k=2,3
+ do l=2,3
+ W1i=W1i+a(k,1)*a(l,1)*Ap_inv(k,l)
+ enddo
+ enddo
+ W1i=a(1,1)-W1i
+ W1(i)=dexp(bsc(i,it)-0.5D0*W1i*zz1*zz1)
+! if (lprint) write(*,'(a,3(1pe15.5)/)')
+! & 'detAp, W1, anormi',detApi,W1i,anormi
+ do k=2,3
+ zk=censc(k,i,it)
+ do l=2,3
+ zk=zk+zz1*Ap_inv(k,l)*a(l,1)
+ enddo
+ z(k,i)=zk
+ enddo
+ detAp(i)=dsqrt(detApi)
+ enddo
+
+ if (lprint) then
+ print *,'W1:',(w1(i),i=1,nlobit)
+ print *,'detAp:',(detAp(i),i=1,nlobit)
+ print *,'Z'
+ do i=1,nlobit
+ print '(i2,3f10.5)',i,(rad2deg*z(j,i),j=2,3)
+ enddo
+ write (iout,*) 'W1:',(w1(i),i=1,nlobit)
+ write (iout,*) 'detAp:',(detAp(i),i=1,nlobit)
+ write (iout,*) 'Z'
+ do i=1,nlobit
+ write (iout,'(i2,3f10.5)') i,(rad2deg*z(j,i),j=2,3)
+ enddo
+ endif
+ if (lcheck) then
+! Writing the distribution just to check the procedure
+ fac=0.0D0
+ dV=deg2rad**2*10.0D0
+ sum=0.0D0
+ sum1=0.0D0
+ do i=1,nlobit
+ fac=fac+W1(i)/detAp(i)
+ enddo
+ fac=1.0D0/(2.0D0*fac*pi)
+!d print *,it,'fac=',fac
+ do ial=90,180,2
+ y(1)=deg2rad*ial
+ do iom=-180,180,5
+ y(2)=deg2rad*iom
+ wart=0.0D0
+ do i=1,nlobit
+ do j=2,3
+ do k=2,3
+ a(j-1,k-1)=gaussc(j,k,i,it)
+ enddo
+ enddo
+ y2=y(2)
+
+ do iii=-1,1
+
+ y(2)=y2+iii*dwapi
+
+ wykl=0.0D0
+ do j=1,2
+ do k=1,2
+ wykl=wykl+a(j,k)*(y(j)-z(j+1,i))*(y(k)-z(k+1,i))
+ enddo
+ enddo
+ wart=wart+W1(i)*dexp(-0.5D0*wykl)
+
+ enddo
+
+ y(2)=y2
+
+ enddo
+! print *,'y',y(1),y(2),' fac=',fac
+ wart=fac*wart
+ write (20,'(2f10.3,1pd15.5)') y(1)*rad2deg,y(2)*rad2deg,wart
+ sum=sum+wart
+ sum1=sum1+1.0D0
+ enddo
+ enddo
+! print *,'it=',it,' sum=',sum*dV,' sum1=',sum1*dV
+ return
+ endif
+
+! Calculate the CM of the system
+!
+ do i=1,nlobit
+ W1(i)=W1(i)/detAp(i)
+ enddo
+ sumW(0)=0.0D0
+ do i=1,nlobit
+ sumW(i)=sumW(i-1)+W1(i)
+ enddo
+ cm(1)=z(2,1)*W1(1)
+ cm(2)=z(3,1)*W1(1)
+ do j=2,nlobit
+ cm(1)=cm(1)+z(2,j)*W1(j)
+ cm(2)=cm(2)+W1(j)*(z(3,1)+pinorm(z(3,j)-z(3,1)))
+ enddo
+ cm(1)=cm(1)/sumW(nlobit)
+ cm(2)=cm(2)/sumW(nlobit)
+ if (cm(1).gt.Big .or. cm(1).lt.-Big .or. &
+ cm(2).gt.Big .or. cm(2).lt.-Big) then
+!d write (iout,'(a)')
+!d & 'Unexpected error in GenSide - CM coordinates too large.'
+!d write (iout,'(i5,2(1pe14.5))') it,cm(1),cm(2)
+!d write (*,'(a)')
+!d & 'Unexpected error in GenSide - CM coordinates too large.'
+!d write (*,'(i5,2(1pe14.5))') it,cm(1),cm(2)
+ fail=.true.
+ return
+ endif
+!d print *,'CM:',cm(1),cm(2)
+!
+! Find the largest search distance from CM
+!
+ radmax=0.0D0
+ do i=1,nlobit
+ do j=2,3
+ do k=2,3
+ a(j-1,k-1)=gaussc(j,k,i,it)
+ enddo
+ enddo
+#ifdef NAG
+ call f02faf('N','U',2,a,3,eig,work,100,ifail)
+#else
+ call djacob(2,3,10000,1.0d-10,a,vec,eig)
+#endif
+#ifdef MPI
+ if (lprint) then
+ print *,'*************** CG Processor',me
+ print *,'CM:',cm(1),cm(2)
+ write (iout,*) '*************** CG Processor',me
+ write (iout,*) 'CM:',cm(1),cm(2)
+ print '(A,8f10.5)','Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
+ write (iout,'(A,8f10.5)') &
+ 'Eigenvalues: ',(1.0/dsqrt(eig(k)),k=1,2)
+ endif
+#endif
+ if (eig(1).lt.eig_limit) then
+ write(iout,'(a)') &
+ 'From Mult_Norm: Eigenvalues of A are too small.'
+ write(*,'(a)') &
+ 'From Mult_Norm: Eigenvalues of A are too small.'
+ fail=.true.
+ return
+ endif
+ radius=0.0D0
+!d print *,'i=',i
+ do j=1,2
+ radius=radius+pinorm(z(j+1,i)-cm(j))**2
+ enddo
+ radius=dsqrt(radius)+3.0D0/dsqrt(eig(1))
+ if (radius.gt.radmax) radmax=radius
+ enddo
+ if (radmax.gt.pi) radmax=pi
+!
+! Determine the boundaries of the search rectangle.
+!
+ if (lprint) then
+ print '(a,4(1pe14.4))','W1: ',(W1(i),i=1,nlob(it) )
+ print '(a,4(1pe14.4))','radmax: ',radmax
+ endif
+ box(1,1)=dmax1(cm(1)-radmax,0.0D0)
+ box(2,1)=dmin1(cm(1)+radmax,pi)
+ box(1,2)=cm(2)-radmax
+ box(2,2)=cm(2)+radmax
+ if (lprint) then
+#ifdef MPI
+ print *,'CG Processor',me,' Array BOX:'
+#else
+ print *,'Array BOX:'
+#endif
+ print '(4(1pe14.4))',((box(k,j),k=1,2),j=1,2)
+ print '(a,4(1pe14.4))','sumW: ',(sumW(i),i=0,nlob(it) )
+#ifdef MPI
+ write (iout,*)'CG Processor',me,' Array BOX:'
+#else
+ write (iout,*)'Array BOX:'
+#endif
+ write(iout,'(4(1pe14.4))') ((box(k,j),k=1,2),j=1,2)
+ write(iout,'(a,4(1pe14.4))')'sumW: ',(sumW(i),i=0,nlob(it) )
+ endif
+ if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then
+#ifdef MPI
+ write (iout,'(a,i4,a,3e15.5)') 'CG Processor:',me,': bad sampling box.',box(1,2),box(2,2),radmax
+ write (*,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
+#else
+! write (iout,'(a)') 'Bad sampling box.'
+#endif
+ fail=.true.
+ return
+ endif
+ which_lobe=ran_number(0.0D0,sumW(nlobit))
+! print '(a,1pe14.4)','which_lobe=',which_lobe
+ do i=1,nlobit
+ if (sumW(i-1).le.which_lobe .and. sumW(i).ge.which_lobe) goto 1
+ enddo
+ 1 ilob=i
+! print *,'ilob=',ilob,' nlob=',nlob(it)
+ do i=2,3
+ cm(i-1)=z(i,ilob)
+ do j=2,3
+ a(i-1,j-1)=gaussc(i,j,ilob,it)
+ enddo
+ enddo
+!d print '(a,i4,a)','CG Processor',me,' Calling MultNorm1.'
+ call mult_norm1(3,2,a,cm,box,y,fail)
+ if (fail) return
+ al=y(1)
+ om=pinorm(y(2))
+!d print *,'al=',al,' om=',om
+!d stop
+ return
+ end subroutine gen_side
+!-----------------------------------------------------------------------------
+ subroutine overlap_sc(scfail)
+!
+! Internal and cartesian coordinates must be consistent as input,
+! and will be up-to-date on return.
+! At the end of this procedure, scfail is true if there are
+! overlapping residues left, or false otherwise (success)
+!
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.FFIELD'
+! include 'COMMON.VAR'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.IOUNITS'
+ logical :: had_overlaps,fail,scfail
+ integer,dimension(nres) :: ioverlap !(maxres)
+ integer :: ioverlap_last,k,maxsi,i,iti,nsi
+ integer :: ires,j
+
+ had_overlaps=.false.
+ call overlap_sc_list(ioverlap,ioverlap_last)
+ if (ioverlap_last.gt.0) then
+ write (iout,*) '#OVERLAPing residues ',ioverlap_last
+ write (iout,'(20i4)') (ioverlap(k),k=1,ioverlap_last)
+ had_overlaps=.true.
+ endif
+
+ maxsi=1000
+ do k=1,1000
+ if (ioverlap_last.eq.0) exit
+
+ do ires=1,ioverlap_last
+ i=ioverlap(ires)
+ iti=iabs(itype(i))
+ if (iti.ne.10) then
+ nsi=0
+ fail=.true.
+ do while (fail.and.nsi.le.maxsi)
+ call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
+ nsi=nsi+1
+ enddo
+ if(fail) goto 999
+ endif
+ enddo
+
+ call chainbuild
+ call overlap_sc_list(ioverlap,ioverlap_last)
+! write (iout,*) 'Overlaping residues ',ioverlap_last,
+! & (ioverlap(j),j=1,ioverlap_last)
+ enddo
+
+ if (k.le.1000.and.ioverlap_last.eq.0) then
+ scfail=.false.
+ if (had_overlaps) then
+ write (iout,*) '#OVERLAPing all corrected after ',k,&
+ ' random generation'
+ endif
+ else
+ scfail=.true.
+ write (iout,*) '#OVERLAPing NOT all corrected ',ioverlap_last
+ write (iout,'(20i4)') (ioverlap(j),j=1,ioverlap_last)
+ endif
+
+ return
+
+ 999 continue
+ write (iout,'(a30,i5,a12,i4)') &
+ '#OVERLAP FAIL in gen_side after',maxsi,&
+ 'iter for RES',i
+ scfail=.true.
+ return
+ end subroutine overlap_sc
+!-----------------------------------------------------------------------------
+ subroutine overlap_sc_list(ioverlap,ioverlap_last)
+ use calc_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.LOCAL'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.FFIELD'
+! include 'COMMON.VAR'
+! include 'COMMON.CALC'
+ logical :: fail
+ integer,dimension(nres) :: ioverlap !(maxres)
+ integer :: ioverlap_last
+!el local variables
+ integer :: ind,iint
+ real(kind=8) :: redfac,sig !rrij,sigsq,
+ integer :: itypi,itypj,itypi1
+ real(kind=8) :: xi,yi,zi,sig0ij,rcomp,rrij,rij_shift
+ data redfac /0.5D0/
+
+ ioverlap_last=0
+! Check for SC-SC overlaps and mark residues
+! print *,'>>overlap_sc nnt=',nnt,' nct=',nct
+ ind=0
+ do i=iatsc_s,iatsc_e
+ itypi=iabs(itype(i))
+ itypi1=iabs(itype(i+1))
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+ dsci_inv=dsc_inv(itypi)
+!
+ do iint=1,nint_gr(i)
+ do j=istart(i,iint),iend(i,iint)
+ ind=ind+1
+ itypj=iabs(itype(j))
+ dscj_inv=dsc_inv(itypj)
+ sig0ij=sigma(itypi,itypj)
+ chi1=chi(itypi,itypj)
+ chi2=chi(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip(itypi)
+ chip2=chip(itypj)
+ chip12=chip1*chip2
+ alf1=alp(itypi)
+ alf2=alp(itypj)
+ alf12=0.5D0*(alf1+alf2)
+ if (j.gt.i+1) then
+ rcomp=sigmaii(itypi,itypj)
+ else
+ rcomp=sigma(itypi,itypj)
+ endif
+! print '(2(a3,2i3),a3,2f10.5)',
+! & ' i=',i,iti,' j=',j,itj,' d=',dist(nres+i,nres+j)
+! & ,rcomp
+ xj=c(1,nres+j)-xi
+ yj=c(2,nres+j)-yi
+ zj=c(3,nres+j)-zi
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+ call sc_angular
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+sig0ij
+
+!t if ( 1.0/rij .lt. redfac*rcomp .or.
+!t & rij_shift.le.0.0D0 ) then
+ if ( rij_shift.le.0.0D0 ) then
+!d write (iout,'(a,i3,a,i3,a,f10.5,a,3f10.5)')
+!d & 'overlap SC-SC: i=',i,' j=',j,
+!d & ' dist=',dist(nres+i,nres+j),' rcomp=',
+!d & rcomp,1.0/rij,rij_shift
+ ioverlap_last=ioverlap_last+1
+ ioverlap(ioverlap_last)=i
+ do k=1,ioverlap_last-1
+ if (ioverlap(k).eq.i) ioverlap_last=ioverlap_last-1
+ enddo
+ ioverlap_last=ioverlap_last+1
+ ioverlap(ioverlap_last)=j
+ do k=1,ioverlap_last-1
+ if (ioverlap(k).eq.j) ioverlap_last=ioverlap_last-1
+ enddo
+ endif
+ enddo
+ enddo
+ enddo
+ return
+ end subroutine overlap_sc_list
+#endif
+!-----------------------------------------------------------------------------
+! energy_p_new_barrier.F
+!-----------------------------------------------------------------------------
+ subroutine sc_angular
+! Calculate eps1,eps2,eps3,sigma, and parts of their derivatives in om1,om2,
+! om12. Called by ebp, egb, and egbv.
+ use calc_data
+! implicit none
+! include 'COMMON.CALC'
+! include 'COMMON.IOUNITS'
+ erij(1)=xj*rij
+ erij(2)=yj*rij
+ erij(3)=zj*rij
+ om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
+ om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
+ om12=dxi*dxj+dyi*dyj+dzi*dzj
+ chiom12=chi12*om12
+! Calculate eps1(om12) and its derivative in om12
+ faceps1=1.0D0-om12*chiom12
+ faceps1_inv=1.0D0/faceps1
+ eps1=dsqrt(faceps1_inv)
+! Following variable is eps1*deps1/dom12
+ eps1_om12=faceps1_inv*chiom12
+! diagnostics only
+! faceps1_inv=om12
+! eps1=om12
+! eps1_om12=1.0d0
+! write (iout,*) "om12",om12," eps1",eps1
+! Calculate sigma(om1,om2,om12) and the derivatives of sigma**2 in om1,om2,
+! and om12.
+ om1om2=om1*om2
+ chiom1=chi1*om1
+ chiom2=chi2*om2
+ facsig=om1*chiom1+om2*chiom2-2.0D0*om1om2*chiom12
+ sigsq=1.0D0-facsig*faceps1_inv
+ sigsq_om1=(chiom1-chiom12*om2)*faceps1_inv
+ sigsq_om2=(chiom2-chiom12*om1)*faceps1_inv
+ sigsq_om12=-chi12*(om1om2*faceps1-om12*facsig)*faceps1_inv**2
+! diagnostics only
+! sigsq=1.0d0
+! sigsq_om1=0.0d0
+! sigsq_om2=0.0d0
+! sigsq_om12=0.0d0
+! write (iout,*) "chiom1",chiom1," chiom2",chiom2," chiom12",chiom12
+! write (iout,*) "faceps1",faceps1," faceps1_inv",faceps1_inv,
+! & " eps1",eps1
+! Calculate eps2 and its derivatives in om1, om2, and om12.
+ chipom1=chip1*om1
+ chipom2=chip2*om2
+ chipom12=chip12*om12
+ facp=1.0D0-om12*chipom12
+ facp_inv=1.0D0/facp
+ facp1=om1*chipom1+om2*chipom2-2.0D0*om1om2*chipom12
+! write (iout,*) "chipom1",chipom1," chipom2",chipom2,
+! & " chipom12",chipom12," facp",facp," facp_inv",facp_inv
+! Following variable is the square root of eps2
+ eps2rt=1.0D0-facp1*facp_inv
+! Following three variables are the derivatives of the square root of eps
+! in om1, om2, and om12.
+ eps2rt_om1=-4.0D0*(chipom1-chipom12*om2)*facp_inv
+ eps2rt_om2=-4.0D0*(chipom2-chipom12*om1)*facp_inv
+ eps2rt_om12=4.0D0*chip12*(om1om2*facp-om12*facp1)*facp_inv**2
+! Evaluate the "asymmetric" factor in the VDW constant, eps3
+ eps3rt=1.0D0-alf1*om1+alf2*om2-alf12*om12
+! write (iout,*) "eps2rt",eps2rt," eps3rt",eps3rt
+! write (iout,*) "eps2rt_om1",eps2rt_om1," eps2rt_om2",eps2rt_om2,
+! & " eps2rt_om12",eps2rt_om12
+! Calculate whole angle-dependent part of epsilon and contributions
+! to its derivatives
+ return
+ end subroutine sc_angular
+!-----------------------------------------------------------------------------
+! initialize_p.F
+!-----------------------------------------------------------------------------
+ subroutine int_bounds(total_ints,lower_bound,upper_bound)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+ include 'mpif.h'
+! include 'COMMON.SETUP'
+ integer :: total_ints,lower_bound,upper_bound,nint
+ integer,dimension(0:nfgtasks) :: int4proc,sint4proc !(0:max_fg_procs)
+ integer :: i,nexcess
+ nint=total_ints/nfgtasks
+ do i=1,nfgtasks
+ int4proc(i-1)=nint
+ enddo
+ nexcess=total_ints-nint*nfgtasks
+ do i=1,nexcess
+ int4proc(nfgtasks-i)=int4proc(nfgtasks-i)+1
+ enddo
+ lower_bound=0
+ do i=0,fg_rank-1
+ lower_bound=lower_bound+int4proc(i)
+ enddo
+ upper_bound=lower_bound+int4proc(fg_rank)
+ lower_bound=lower_bound+1
+ return
+ end subroutine int_bounds
+!-----------------------------------------------------------------------------
+ subroutine int_bounds1(total_ints,lower_bound,upper_bound)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+ include 'mpif.h'
+! include 'COMMON.SETUP'
+ integer :: total_ints,lower_bound,upper_bound,nint
+ integer :: nexcess,i
+ integer,dimension(0:nfgtasks) :: int4proc,sint4proc !(0:max_fg_procs)
+ nint=total_ints/nfgtasks1
+ do i=1,nfgtasks1
+ int4proc(i-1)=nint
+ enddo
+ nexcess=total_ints-nint*nfgtasks1
+ do i=1,nexcess
+ int4proc(nfgtasks1-i)=int4proc(nfgtasks1-i)+1
+ enddo
+ lower_bound=0
+ do i=0,fg_rank1-1
+ lower_bound=lower_bound+int4proc(i)
+ enddo
+ upper_bound=lower_bound+int4proc(fg_rank1)
+ lower_bound=lower_bound+1
+ return
+ end subroutine int_bounds1
+!-----------------------------------------------------------------------------
+! intcartderiv.F
+!-----------------------------------------------------------------------------
+ subroutine chainbuild_cart
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+ use control_data
+#ifdef MPI
+ include 'mpif.h'
+#endif
+! include 'COMMON.SETUP'
+! include 'COMMON.CHAIN'
+! include 'COMMON.LOCAL'
+! include 'COMMON.TIME1'
+! include 'COMMON.IOUNITS'
+ integer :: j,i,ierror,ierr
+ real(kind=8) :: time00,time01
+#ifdef MPI
+ if (nfgtasks.gt.1) then
+! write (iout,*) "BCAST in chainbuild_cart"
+! call flush(iout)
+! Broadcast the order to build the chain and compute internal coordinates
+! to the slaves. The slaves receive the order in ERGASTULUM.
+ time00=MPI_Wtime()
+! write (iout,*) "CHAINBUILD_CART: DC before BCAST"
+! do i=0,nres
+! write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
+! & (dc(j,i+nres),j=1,3)
+! enddo
+ if (fg_rank.eq.0) &
+ call MPI_Bcast(7,1,MPI_INTEGER,king,FG_COMM,IERROR)
+ time_bcast7=time_bcast7+MPI_Wtime()-time00
+ time01=MPI_Wtime()
+ call MPI_Bcast(dc(1,0),6*(nres+1),MPI_DOUBLE_PRECISION,&
+ king,FG_COMM,IERR)
+! write (iout,*) "CHAINBUILD_CART: DC after BCAST"
+! do i=0,nres
+! write (iout,'(i3,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
+! & (dc(j,i+nres),j=1,3)
+! enddo
+! write (iout,*) "End BCAST in chainbuild_cart"
+! call flush(iout)
+ time_bcast=time_bcast+MPI_Wtime()-time00
+ time_bcastc=time_bcastc+MPI_Wtime()-time01
+ endif
+#endif
+ do j=1,3
+ c(j,1)=dc(j,0)
+ enddo
+ do i=2,nres
+ do j=1,3
+ c(j,i)=c(j,i-1)+dc(j,i-1)
+ enddo
+ enddo
+ do i=1,nres
+ do j=1,3
+ c(j,i+nres)=c(j,i)+dc(j,i+nres)
+ enddo
+ enddo
+! write (iout,*) "CHAINBUILD_CART"
+! call cartprint
+ call int_from_cart1(.false.)
+ return
+ end subroutine chainbuild_cart
+!-----------------------------------------------------------------------------
+! intcor.f
+!-----------------------------------------------------------------------------
+ real(kind=8) function alpha(i1,i2,i3)
+!
+! Calculates the planar angle between atoms (i1), (i2), and (i3).
+!
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.CHAIN'
+!el local variables
+ integer :: i1,i2,i3
+ real(kind=8) :: x12,x23,y12,y23,z12,z23,vnorm,wnorm,scalar
+ x12=c(1,i1)-c(1,i2)
+ x23=c(1,i3)-c(1,i2)
+ y12=c(2,i1)-c(2,i2)
+ y23=c(2,i3)-c(2,i2)
+ z12=c(3,i1)-c(3,i2)
+ z23=c(3,i3)-c(3,i2)
+ vnorm=dsqrt(x12*x12+y12*y12+z12*z12)
+ wnorm=dsqrt(x23*x23+y23*y23+z23*z23)
+ scalar=(x12*x23+y12*y23+z12*z23)/(vnorm*wnorm)
+ alpha=arcos(scalar)
+ return
+ end function alpha
+!-----------------------------------------------------------------------------
+ real(kind=8) function beta(i1,i2,i3,i4)
+!
+! Calculates the dihedral angle between atoms (i1), (i2), (i3) and (i4)
+!
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.CHAIN'
+!el local variables
+ integer :: i1,i2,i3,i4
+ real(kind=8) :: x12,x23,x34,y12,y23,y34,z12,z23,z34
+ real(kind=8) :: wx,wy,wz,wnorm,vx,vy,vz,vnorm,scalar,angle
+ real(kind=8) :: tx,ty,tz
+ x12=c(1,i1)-c(1,i2)
+ x23=c(1,i3)-c(1,i2)
+ x34=c(1,i4)-c(1,i3)
+ y12=c(2,i1)-c(2,i2)
+ y23=c(2,i3)-c(2,i2)
+ y34=c(2,i4)-c(2,i3)
+ z12=c(3,i1)-c(3,i2)
+ z23=c(3,i3)-c(3,i2)
+ z34=c(3,i4)-c(3,i3)
+!d print '(2i3,3f10.5)',i1,i2,x12,y12,z12
+!d print '(2i3,3f10.5)',i2,i3,x23,y23,z23
+!d print '(2i3,3f10.5)',i3,i4,x34,y34,z34
+ wx=-y23*z34+y34*z23
+ wy=x23*z34-z23*x34
+ wz=-x23*y34+y23*x34
+ wnorm=dsqrt(wx*wx+wy*wy+wz*wz)
+ vx=y12*z23-z12*y23
+ vy=-x12*z23+z12*x23
+ vz=x12*y23-y12*x23
+ vnorm=dsqrt(vx*vx+vy*vy+vz*vz)
+ if (vnorm.gt.1.0D-13 .and. wnorm.gt.1.0D-13) then
+ scalar=(vx*wx+vy*wy+vz*wz)/(vnorm*wnorm)
+ if (dabs(scalar).gt.1.0D0) &
+ scalar=0.99999999999999D0*scalar/dabs(scalar)
+ angle=dacos(scalar)
+!d print '(2i4,10f7.3)',i2,i3,vx,vy,vz,wx,wy,wz,vnorm,wnorm,
+!d &scalar,angle
+ else
+ angle=pi
+ endif
+! if (angle.le.0.0D0) angle=pi+angle
+ tx=vy*wz-vz*wy
+ ty=-vx*wz+vz*wx
+ tz=vx*wy-vy*wx
+ scalar=tx*x23+ty*y23+tz*z23
+ if (scalar.lt.0.0D0) angle=-angle
+ beta=angle
+ return
+ end function beta
+!-----------------------------------------------------------------------------
+ real(kind=8) function dist(i1,i2)
+!
+! Calculates the distance between atoms (i1) and (i2).
+!
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.CHAIN'
+!el local variables
+ integer :: i1,i2
+ real(kind=8) :: x12,y12,z12
+ x12=c(1,i1)-c(1,i2)
+ y12=c(2,i1)-c(2,i2)
+ z12=c(3,i1)-c(3,i2)
+ dist=dsqrt(x12*x12+y12*y12+z12*z12)
+ return
+ end function dist
+#if !defined(WHAM_RUN) && !defined(CLUSTER)
+!-----------------------------------------------------------------------------
+! local_move.f
+!-----------------------------------------------------------------------------
+ subroutine local_move_init(debug)
+!rc implicit none
+
+! Includes
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS' ! Needed by COMMON.LOCAL
+! include 'COMMON.GEO' ! For pi, deg2rad
+! include 'COMMON.LOCAL' ! For vbl
+! include 'COMMON.LOCMOVE'
+
+! INPUT arguments
+ logical :: debug
+
+
+! Determine wheter to do some debugging output
+ locmove_output=debug
+
+! Set the init_called flag to 1
+ init_called=1
+
+! The following are never changed
+ min_theta=60.D0*deg2rad ! (0,PI)
+ max_theta=175.D0*deg2rad ! (0,PI)
+ dmin2=vbl*vbl*2.*(1.-cos(min_theta))
+ dmax2=vbl*vbl*2.*(1.-cos(max_theta))
+ flag=1.0D300
+ small=1.0D-5
+ small2=0.5*small*small
+
+! Not really necessary...
+ a_n=0
+ b_n=0
+ res_n=0
+
+ return
+ end subroutine local_move_init
+!-----------------------------------------------------------------------------
+ subroutine local_move(n_start, n_end, PHImin, PHImax)
+! Perform a local move between residues m and n (inclusive)
+! PHImin and PHImax [0,PI] determine the size of the move
+! Works on whatever structure is in the variables theta and phi,
+! sidechain variables are left untouched
+! The final structure is NOT minimized, but both the cartesian
+! variables c and the angles are up-to-date at the end (no further
+! chainbuild is required)
+!rc implicit none
+ use random,only:ran_number
+! Includes
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.CHAIN'
+! include 'COMMON.VAR'
+! include 'COMMON.MINIM'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.LOCMOVE'
+
+! External functions
+!EL integer move_res
+!EL external move_res
+!EL double precision ran_number
+!EL external ran_number
+
+! INPUT arguments
+ integer :: n_start, n_end ! First and last residues to move
+ real(kind=8) :: PHImin, PHImax ! min/max angles [0,PI]
+
+! Local variables
+ integer :: i,j
+ real(kind=8) :: min,max
+ integer :: iretcode
+
+
+! Check if local_move_init was called. This assumes that it
+! would not be 1 if not explicitely initialized
+ if (init_called.ne.1) then
+ write(6,*)' *** local_move_init not called!!!'
+ stop
+ endif
+
+! Quick check for crazy range
+ if (n_start.gt.n_end .or. n_start.lt.1 .or. n_end.gt.nres) then
+ write(6,'(a,i3,a,i3)') &
+ ' *** Cannot make local move between n_start = ',&
+ n_start,' and n_end = ',n_end
+ return
+ endif
+
+! Take care of end residues first...
+ if (n_start.eq.1) then
+! Move residue 1 (completely random)
+ theta(3)=ran_number(min_theta,max_theta)
+ phi(4)=ran_number(-PI,PI)
+ i=2
+ else
+ i=n_start
+ endif
+ if (n_end.eq.nres) then
+! Move residue nres (completely random)
+ theta(nres)=ran_number(min_theta,max_theta)
+ phi(nres)=ran_number(-PI,PI)
+ j=nres-1
+ else
+ j=n_end
+ endif
+
+! ...then go through all other residues one by one
+! Start from the two extremes and converge
+ call chainbuild
+ do while (i.le.j)
+ min=PHImin
+ max=PHImax
+!$$$c Move the first two residues by less than the others
+!$$$ if (i-n_start.lt.3) then
+!$$$ if (i-n_start.eq.0) then
+!$$$ min=0.4*PHImin
+!$$$ max=0.4*PHImax
+!$$$ else if (i-n_start.eq.1) then
+!$$$ min=0.8*PHImin
+!$$$ max=0.8*PHImax
+!$$$ else if (i-n_start.eq.2) then
+!$$$ min=PHImin
+!$$$ max=PHImax
+!$$$ endif
+!$$$ endif
+
+! The actual move, on residue i
+ iretcode=move_res(min,max,i) ! Discard iretcode
+ i=i+1
+
+ if (i.le.j) then
+ min=PHImin
+ max=PHImax
+!$$$c Move the last two residues by less than the others
+!$$$ if (n_end-j.lt.3) then
+!$$$ if (n_end-j.eq.0) then
+!$$$ min=0.4*PHImin
+!$$$ max=0.4*PHImax
+!$$$ else if (n_end-j.eq.1) then
+!$$$ min=0.8*PHImin
+!$$$ max=0.8*PHImax
+!$$$ else if (n_end-j.eq.2) then
+!$$$ min=PHImin
+!$$$ max=PHImax
+!$$$ endif
+!$$$ endif
+
+! The actual move, on residue j
+ iretcode=move_res(min,max,j) ! Discard iretcode
+ j=j-1
+ endif
+ enddo
+
+ call int_from_cart(.false.,.false.)
+
+ return
+ end subroutine local_move
+!-----------------------------------------------------------------------------
+ subroutine output_tabs
+! Prints out the contents of a_..., b_..., res_...
+! implicit none
+
+! Includes
+! include 'COMMON.GEO'
+! include 'COMMON.LOCMOVE'
+
+! Local variables
+ integer :: i,j
+
+ write(6,*)'a_...'
+ write(6,'(8f7.1)')(a_ang(i)*rad2deg,i=0,a_n-1)
+ write(6,'(8(2x,3l1,2x))')((a_tab(i,j),i=0,2),j=0,a_n-1)
+
+ write(6,*)'b_...'
+ write(6,'(4f7.1)')(b_ang(i)*rad2deg,i=0,b_n-1)
+ write(6,'(4(2x,3l1,2x))')((b_tab(i,j),i=0,2),j=0,b_n-1)
+
+ write(6,*)'res_...'
+ write(6,'(12f7.1)')(res_ang(i)*rad2deg,i=0,res_n-1)
+ write(6,'(12(2x,3l1,2x))')((res_tab(0,i,j),i=0,2),j=0,res_n-1)
+ write(6,'(12(2x,3l1,2x))')((res_tab(1,i,j),i=0,2),j=0,res_n-1)
+ write(6,'(12(2x,3l1,2x))')((res_tab(2,i,j),i=0,2),j=0,res_n-1)
+
+ return
+ end subroutine output_tabs
+!-----------------------------------------------------------------------------
+ subroutine angles2tab(PHImin,PHImax,n,ang,tab)
+! Only uses angles if [0,PI] (but PHImin cannot be 0.,
+! and PHImax cannot be PI)
+! implicit none
+
+! Includes
+! include 'COMMON.GEO'
+
+! INPUT arguments
+ real(kind=8) :: PHImin,PHImax
+
+! OUTPUT arguments
+ integer :: n
+ real(kind=8),dimension(0:3) :: ang
+ logical,dimension(0:2,0:3) :: tab
+
+
+ if (PHImin .eq. PHImax) then
+! Special case with two 010's
+ n = 2;
+ ang(0) = -PHImin;
+ ang(1) = PHImin;
+ tab(0,0) = .false.
+ tab(2,0) = .false.
+ tab(0,1) = .false.
+ tab(2,1) = .false.
+ tab(1,0) = .true.
+ tab(1,1) = .true.
+ else if (PHImin .eq. PI) then
+! Special case with one 010
+ n = 1
+ ang(0) = PI
+ tab(0,0) = .false.
+ tab(2,0) = .false.
+ tab(1,0) = .true.
+ else if (PHImax .eq. 0.) then
+! Special case with one 010
+ n = 1
+ ang(0) = 0.
+ tab(0,0) = .false.
+ tab(2,0) = .false.
+ tab(1,0) = .true.
+ else
+! Standard cases
+ n = 0
+ if (PHImin .gt. 0.) then
+! Start of range (011)
+ ang(n) = PHImin
+ tab(0,n) = .false.
+ tab(1,n) = .true.
+ tab(2,n) = .true.
+! End of range (110)
+ ang(n+1) = -PHImin
+ tab(0,n+1) = .true.
+ tab(1,n+1) = .true.
+ tab(2,n+1) = .false.
+ n = n+2
+ endif
+ if (PHImax .lt. PI) then
+! Start of range (011)
+ ang(n) = -PHImax
+ tab(0,n) = .false.
+ tab(1,n) = .true.
+ tab(2,n) = .true.
+! End of range (110)
+ ang(n+1) = PHImax
+ tab(0,n+1) = .true.
+ tab(1,n+1) = .true.
+ tab(2,n+1) = .false.
+ n = n+2
+ endif
+ endif
+
+ return
+ end subroutine angles2tab
+!-----------------------------------------------------------------------------
+ subroutine minmax_angles(x,y,z,r,n,ang,tab)
+! When solutions do not exist, assume all angles
+! are acceptable - i.e., initial geometry must be correct
+! implicit none
+
+! Includes
+! include 'COMMON.GEO'
+! include 'COMMON.LOCMOVE'
+
+! Input arguments
+ real(kind=8) :: x,y,z,r
+
+! Output arguments
+ integer :: n
+ real(kind=8),dimension(0:3) :: ang
+ logical,dimension(0:2,0:3) :: tab
+
+! Local variables
+ real(kind=8) :: num, denom, phi
+ real(kind=8) :: Kmin, Kmax
+ integer :: i
+
+
+ num = x*x + y*y + z*z
+ denom = x*x + y*y
+ n = 0
+ if (denom .gt. 0.) then
+ phi = atan2(y,x)
+ denom = 2.*r*sqrt(denom)
+ num = num+r*r
+ Kmin = (num - dmin2)/denom
+ Kmax = (num - dmax2)/denom
+
+! Allowed values of K (else all angles are acceptable)
+! -1 <= Kmin < 1
+! -1 < Kmax <= 1
+ if (Kmin .gt. 1. .or. abs(Kmin-1.) .lt. small2) then
+ Kmin = -flag
+ else if (Kmin .lt. -1. .or. abs(Kmin+1.) .lt. small2) then
+ Kmin = PI
+ else
+ Kmin = acos(Kmin)
+ endif
+
+ if (Kmax .lt. -1. .or. abs(Kmax+1.) .lt. small2) then
+ Kmax = flag
+ else if (Kmax .gt. 1. .or. abs(Kmax-1.) .lt. small2) then
+ Kmax = 0.
+ else
+ Kmax = acos(Kmax)
+ endif
+
+ if (Kmax .lt. Kmin) Kmax = Kmin
+
+ call angles2tab(Kmin, Kmax, n, ang, tab)
+
+! Add phi and check that angles are within range (-PI,PI]
+ do i=0,n-1
+ ang(i) = ang(i)+phi
+ if (ang(i) .le. -PI) then
+ ang(i) = ang(i)+2.*PI
+ else if (ang(i) .gt. PI) then
+ ang(i) = ang(i)-2.*PI
+ endif
+ enddo
+ endif
+
+ return
+ end subroutine minmax_angles
+!-----------------------------------------------------------------------------
+ subroutine construct_tab
+! Take a_... and b_... values and produces the results res_...
+! x_ang are assumed to be all different (diff > small)
+! x_tab(1,i) must be 1 for all i (i.e., all x_ang are acceptable)
+! implicit none
+
+! Includes
+! include 'COMMON.LOCMOVE'
+
+! Local variables
+ integer :: n_max,i,j,index
+ logical :: done
+ real(kind=8) :: phi
+
+
+ n_max = a_n + b_n
+ if (n_max .eq. 0) then
+ res_n = 0
+ return
+ endif
+
+ do i=0,n_max-1
+ do j=0,1
+ res_tab(j,0,i) = .true.
+ res_tab(j,2,i) = .true.
+ res_tab(j,1,i) = .false.
+ enddo
+ enddo
+
+ index = 0
+ phi = -flag
+ done = .false.
+ do while (.not.done)
+ res_ang(index) = flag
+
+! Check a first...
+ do i=0,a_n-1
+ if ((a_ang(i)-phi).gt.small .and. &
+ a_ang(i) .lt. res_ang(index)) then
+! Found a lower angle
+ res_ang(index) = a_ang(i)
+! Copy the values from a_tab into res_tab(0,,)
+ res_tab(0,0,index) = a_tab(0,i)
+ res_tab(0,1,index) = a_tab(1,i)
+ res_tab(0,2,index) = a_tab(2,i)
+! Set default values for res_tab(1,,)
+ res_tab(1,0,index) = .true.
+ res_tab(1,1,index) = .false.
+ res_tab(1,2,index) = .true.
+ else if (abs(a_ang(i)-res_ang(index)).lt.small) then
+! Found an equal angle (can only be equal to a b_ang)
+ res_tab(0,0,index) = a_tab(0,i)
+ res_tab(0,1,index) = a_tab(1,i)
+ res_tab(0,2,index) = a_tab(2,i)
+ endif
+ enddo
+! ...then check b
+ do i=0,b_n-1
+ if ((b_ang(i)-phi).gt.small .and. &
+ b_ang(i) .lt. res_ang(index)) then
+! Found a lower angle
+ res_ang(index) = b_ang(i)
+! Copy the values from b_tab into res_tab(1,,)
+ res_tab(1,0,index) = b_tab(0,i)
+ res_tab(1,1,index) = b_tab(1,i)
+ res_tab(1,2,index) = b_tab(2,i)
+! Set default values for res_tab(0,,)
+ res_tab(0,0,index) = .true.
+ res_tab(0,1,index) = .false.
+ res_tab(0,2,index) = .true.
+ else if (abs(b_ang(i)-res_ang(index)).lt.small) then
+! Found an equal angle (can only be equal to an a_ang)
+ res_tab(1,0,index) = b_tab(0,i)
+ res_tab(1,1,index) = b_tab(1,i)
+ res_tab(1,2,index) = b_tab(2,i)
+ endif
+ enddo
+
+ if (res_ang(index) .eq. flag) then
+ res_n = index
+ done = .true.
+ else if (index .eq. n_max-1) then
+ res_n = n_max
+ done = .true.
+ else
+ phi = res_ang(index) ! Store previous angle
+ index = index+1
+ endif
+ enddo
+
+! Fill the gaps
+! First a...
+ index = 0
+ if (a_n .gt. 0) then
+ do while (.not.res_tab(0,1,index))
+ index=index+1
+ enddo
+ done = res_tab(0,2,index)
+ do i=index+1,res_n-1
+ if (res_tab(0,1,i)) then
+ done = res_tab(0,2,i)
+ else
+ res_tab(0,0,i) = done
+ res_tab(0,1,i) = done
+ res_tab(0,2,i) = done
+ endif
+ enddo
+ done = res_tab(0,0,index)
+ do i=index-1,0,-1
+ if (res_tab(0,1,i)) then
+ done = res_tab(0,0,i)
+ else
+ res_tab(0,0,i) = done
+ res_tab(0,1,i) = done
+ res_tab(0,2,i) = done
+ endif
+ enddo
+ else
+ do i=0,res_n-1
+ res_tab(0,0,i) = .true.
+ res_tab(0,1,i) = .true.
+ res_tab(0,2,i) = .true.
+ enddo
+ endif
+! ...then b
+ index = 0
+ if (b_n .gt. 0) then
+ do while (.not.res_tab(1,1,index))
+ index=index+1
+ enddo
+ done = res_tab(1,2,index)
+ do i=index+1,res_n-1
+ if (res_tab(1,1,i)) then
+ done = res_tab(1,2,i)
+ else
+ res_tab(1,0,i) = done
+ res_tab(1,1,i) = done
+ res_tab(1,2,i) = done
+ endif
+ enddo
+ done = res_tab(1,0,index)
+ do i=index-1,0,-1
+ if (res_tab(1,1,i)) then
+ done = res_tab(1,0,i)
+ else
+ res_tab(1,0,i) = done
+ res_tab(1,1,i) = done
+ res_tab(1,2,i) = done
+ endif
+ enddo
+ else
+ do i=0,res_n-1
+ res_tab(1,0,i) = .true.
+ res_tab(1,1,i) = .true.
+ res_tab(1,2,i) = .true.
+ enddo
+ endif
+
+! Finally fill the last row with AND operation
+ do i=0,res_n-1
+ do j=0,2
+ res_tab(2,j,i) = (res_tab(0,j,i) .and. res_tab(1,j,i))
+ enddo
+ enddo
+
+ return
+ end subroutine construct_tab
+!-----------------------------------------------------------------------------
+ subroutine construct_ranges(phi_n,phi_start,phi_end)
+! Given the data in res_..., construct a table of
+! min/max allowed angles
+! implicit none
+
+! Includes
+! include 'COMMON.GEO'
+! include 'COMMON.LOCMOVE'
+
+! Output arguments
+ integer :: phi_n
+ real(kind=8),dimension(0:11) :: phi_start,phi_end
+
+! Local variables
+ logical :: done
+ integer :: index
+
+
+ if (res_n .eq. 0) then
+! Any move is allowed
+ phi_n = 1
+ phi_start(0) = -PI
+ phi_end(0) = PI
+ else
+ phi_n = 0
+ index = 0
+ done = .false.
+ do while (.not.done)
+! Find start of range (01x)
+ done = .false.
+ do while (.not.done)
+ if (res_tab(2,0,index).or.(.not.res_tab(2,1,index))) then
+ index=index+1
+ else
+ done = .true.
+ phi_start(phi_n) = res_ang(index)
+ endif
+ if (index .eq. res_n) done = .true.
+ enddo
+! If a start was found (index < res_n), find the end of range (x10)
+! It may not be found without wrapping around
+ if (index .lt. res_n) then
+ done = .false.
+ do while (.not.done)
+ if ((.not.res_tab(2,1,index)).or.res_tab(2,2,index)) then
+ index=index+1
+ else
+ done = .true.
+ endif
+ if (index .eq. res_n) done = .true.
+ enddo
+ if (index .lt. res_n) then
+! Found the end of the range
+ phi_end(phi_n) = res_ang(index)
+ phi_n=phi_n+1
+ index=index+1
+ if (index .eq. res_n) then
+ done = .true.
+ else
+ done = .false.
+ endif
+ else
+! Need to wrap around
+ done = .true.
+ phi_end(phi_n) = flag
+ endif
+ endif
+ enddo
+! Take care of the last one if need to wrap around
+ if (phi_end(phi_n) .eq. flag) then
+ index = 0
+ do while ((.not.res_tab(2,1,index)).or.res_tab(2,2,index))
+ index=index+1
+ enddo
+ phi_end(phi_n) = res_ang(index) + 2.*PI
+ phi_n=phi_n+1
+ endif
+ endif
+
+ return
+ end subroutine construct_ranges
+!-----------------------------------------------------------------------------
+ subroutine fix_no_moves(phi)
+! implicit none
+
+! Includes
+! include 'COMMON.GEO'
+! include 'COMMON.LOCMOVE'
+
+! Output arguments
+ real(kind=8) :: phi
+
+! Local variables
+ integer :: index
+ real(kind=8) :: diff,temp
+
+
+! Look for first 01x in gammas (there MUST be at least one)
+ diff = flag
+ index = 0
+ do while (res_tab(1,0,index) .or. (.not.res_tab(1,1,index)))
+ index=index+1
+ enddo
+ if (res_ang(index) .le. 0.D0) then ! Make sure it's from PHImax
+! Try to increase PHImax
+ if (index .gt. 0) then
+ phi = res_ang(index-1)
+ diff = abs(res_ang(index) - res_ang(index-1))
+ endif
+! Look for last (corresponding) x10
+ index = res_n - 1
+ do while ((.not.res_tab(1,1,index)) .or. res_tab(1,2,index))
+ index=index-1
+ enddo
+ if (index .lt. res_n-1) then
+ temp = abs(res_ang(index) - res_ang(index+1))
+ if (temp .lt. diff) then
+ phi = res_ang(index+1)
+ diff = temp
+ endif
+ endif
+ endif
+
+! If increasing PHImax didn't work, decreasing PHImin
+! will (with one exception)
+! Look for first x10 (there MUST be at least one)
+ index = 0
+ do while ((.not.res_tab(1,1,index)) .or. res_tab(1,2,index))
+ index=index+1
+ enddo
+ if (res_ang(index) .lt. 0.D0) then ! Make sure it's from PHImin
+! Try to decrease PHImin
+ if (index .lt. res_n-1) then
+ temp = abs(res_ang(index) - res_ang(index+1))
+ if (res_ang(index+1) .le. 0.D0 .and. temp .lt. diff) then
+ phi = res_ang(index+1)
+ diff = temp
+ endif
+ endif
+! Look for last (corresponding) 01x
+ index = res_n - 1
+ do while (res_tab(1,0,index) .or. (.not.res_tab(1,1,index)))
+ index=index-1
+ enddo
+ if (index .gt. 0) then
+ temp = abs(res_ang(index) - res_ang(index-1))
+ if (res_ang(index-1) .ge. 0.D0 .and. temp .lt. diff) then
+ phi = res_ang(index-1)
+ diff = temp
+ endif
+ endif
+ endif
+
+! If it still didn't work, it must be PHImax == 0. or PHImin == PI
+ if (diff .eq. flag) then
+ index = 0
+ if (res_tab(index,1,0) .or. (.not.res_tab(index,1,1)) .or. &
+ res_tab(index,1,2)) index = res_n - 1
+! This MUST work at this point
+ if (index .eq. 0) then
+ phi = res_ang(1)
+ else
+ phi = res_ang(index - 1)
+ endif
+ endif
+
+ return
+ end subroutine fix_no_moves
+!-----------------------------------------------------------------------------
+ integer function move_res(PHImin,PHImax,i_move)
+! Moves residue i_move (in array c), leaving everything else fixed
+! Starting geometry is not checked, it should be correct!
+! R(,i_move) is the only residue that will move, but must have
+! 1 < i_move < nres (i.e., cannot move ends)
+! Whether any output is done is controlled by locmove_output
+!rc implicit none
+ use random,only:ran_number
+! Includes
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! include 'COMMON.LOCMOVE'
+
+! External functions
+!EL double precision ran_number
+!EL external ran_number
+
+! Input arguments
+ real(kind=8) :: PHImin,PHImax
+ integer :: i_move
+
+! RETURN VALUES:
+! 0: move successfull
+! 1: Dmin or Dmax had to be modified
+! 2: move failed - check your input geometry
+
+
+! Local variables
+ real(kind=8),dimension(0:2) :: X,Y,Z,Orig
+ real(kind=8),dimension(0:2) :: P
+ logical :: no_moves,done
+ integer :: index,i,j
+ real(kind=8) :: phi,temp,radius
+ real(kind=8),dimension(0:11) :: phi_start,phi_end
+ integer :: phi_n
+
+! Set up the coordinate system
+ do i=0,2
+ Orig(i)=0.5*(c(i+1,i_move-1)+c(i+1,i_move+1)) ! Position of origin
+ enddo
+
+ do i=0,2
+ Z(i)=c(i+1,i_move+1)-c(i+1,i_move-1)
+ enddo
+ temp=sqrt(Z(0)*Z(0)+Z(1)*Z(1)+Z(2)*Z(2))
+ do i=0,2
+ Z(i)=Z(i)/temp
+ enddo
+
+ do i=0,2
+ X(i)=c(i+1,i_move)-Orig(i)
+ enddo
+! radius is the radius of the circle on which c(,i_move) can move
+ radius=sqrt(X(0)*X(0)+X(1)*X(1)+X(2)*X(2))
+ do i=0,2
+ X(i)=X(i)/radius
+ enddo
+
+ Y(0)=Z(1)*X(2)-X(1)*Z(2)
+ Y(1)=X(0)*Z(2)-Z(0)*X(2)
+ Y(2)=Z(0)*X(1)-X(0)*Z(1)
+
+! Calculate min, max angles coming from dmin, dmax to c(,i_move-2)
+ if (i_move.gt.2) then
+ do i=0,2
+ P(i)=c(i+1,i_move-2)-Orig(i)
+ enddo
+ call minmax_angles(P(0)*X(0)+P(1)*X(1)+P(2)*X(2),&
+ P(0)*Y(0)+P(1)*Y(1)+P(2)*Y(2),&
+ P(0)*Z(0)+P(1)*Z(1)+P(2)*Z(2),&
+ radius,a_n,a_ang,a_tab)
+ else
+ a_n=0
+ endif
+
+! Calculate min, max angles coming from dmin, dmax to c(,i_move+2)
+ if (i_move.lt.nres-2) then
+ do i=0,2
+ P(i)=c(i+1,i_move+2)-Orig(i)
+ enddo
+ call minmax_angles(P(0)*X(0)+P(1)*X(1)+P(2)*X(2),&
+ P(0)*Y(0)+P(1)*Y(1)+P(2)*Y(2),&
+ P(0)*Z(0)+P(1)*Z(1)+P(2)*Z(2),&
+ radius,b_n,b_ang,b_tab)
+ else
+ b_n=0
+ endif
+
+! Construct the resulting table for alpha and beta
+ call construct_tab()
+
+ if (locmove_output) then
+ print *,'ALPHAS & BETAS TABLE'
+ call output_tabs()
+ endif
+
+! Check that there is at least one possible move
+ no_moves = .true.
+ if (res_n .eq. 0) then
+ no_moves = .false.
+ else
+ index = 0
+ do while ((index .lt. res_n) .and. no_moves)
+ if (res_tab(2,1,index)) no_moves = .false.
+ index=index+1
+ enddo
+ endif
+ if (no_moves) then
+ if (locmove_output) print *,' *** Cannot move anywhere'
+ move_res=2
+ return
+ endif
+
+! Transfer res_... into a_...
+ a_n = 0
+ do i=0,res_n-1
+ if ( (res_tab(2,0,i).neqv.res_tab(2,1,i)) .or. &
+ (res_tab(2,0,i).neqv.res_tab(2,2,i)) ) then
+ a_ang(a_n) = res_ang(i)
+ do j=0,2
+ a_tab(j,a_n) = res_tab(2,j,i)
+ enddo
+ a_n=a_n+1
+ endif
+ enddo
+
+! Check that the PHI's are within [0,PI]
+ if (PHImin .lt. 0. .or. abs(PHImin) .lt. small) PHImin = -flag
+ if (PHImin .gt. PI .or. abs(PHImin-PI) .lt. small) PHImin = PI
+ if (PHImax .gt. PI .or. abs(PHImax-PI) .lt. small) PHImax = flag
+ if (PHImax .lt. 0. .or. abs(PHImax) .lt. small) PHImax = 0.
+ if (PHImax .lt. PHImin) PHImax = PHImin
+! Calculate min and max angles coming from PHImin and PHImax,
+! and put them in b_...
+ call angles2tab(PHImin, PHImax, b_n, b_ang, b_tab)
+! Construct the final table
+ call construct_tab()
+
+ if (locmove_output) then
+ print *,'FINAL TABLE'
+ call output_tabs()
+ endif
+
+! Check that there is at least one possible move
+ no_moves = .true.
+ if (res_n .eq. 0) then
+ no_moves = .false.
+ else
+ index = 0
+ do while ((index .lt. res_n) .and. no_moves)
+ if (res_tab(2,1,index)) no_moves = .false.
+ index=index+1
+ enddo
+ endif
+
+ if (no_moves) then
+! Take care of the case where no solution exists...
+ call fix_no_moves(phi)
+ if (locmove_output) then
+ print *,' *** Had to modify PHImin or PHImax'
+ print *,'phi: ',phi*rad2deg
+ endif
+ move_res=1
+ else
+! ...or calculate the solution
+! Construct phi_start/phi_end arrays
+ call construct_ranges(phi_n, phi_start, phi_end)
+! Choose random angle phi in allowed range(s)
+ temp = 0.
+ do i=0,phi_n-1
+ temp = temp + phi_end(i) - phi_start(i)
+ enddo
+ phi = ran_number(phi_start(0),phi_start(0)+temp)
+ index = 0
+ done = .false.
+ do while (.not.done)
+ if (phi .lt. phi_end(index)) then
+ done = .true.
+ else
+ index=index+1
+ endif
+ if (index .eq. phi_n) then
+ done = .true.
+ else if (.not.done) then
+ phi = phi + phi_start(index) - phi_end(index-1)
+ endif
+ enddo
+ if (index.eq.phi_n) phi=phi_end(phi_n-1) ! Fix numerical errors
+ if (phi .gt. PI) phi = phi-2.*PI
+
+ if (locmove_output) then
+ print *,'ALLOWED RANGE(S)'
+ do i=0,phi_n-1
+ print *,phi_start(i)*rad2deg,phi_end(i)*rad2deg
+ enddo
+ print *,'phi: ',phi*rad2deg
+ endif
+ move_res=0
+ endif
+
+! Re-use radius as temp variable
+ temp=radius*cos(phi)
+ radius=radius*sin(phi)
+ do i=0,2
+ c(i+1,i_move)=Orig(i)+temp*X(i)+radius*Y(i)
+ enddo
+
+ return
+ end function move_res
+!-----------------------------------------------------------------------------
+ subroutine loc_test
+!rc implicit none
+
+! Includes
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.LOCAL'
+! include 'COMMON.LOCMOVE'
+
+! External functions
+!EL integer move_res
+!EL external move_res
+
+! Local variables
+ integer :: i,j,imov
+ integer :: phi_n
+ real(kind=8),dimension(0:11) :: phi_start,phi_end
+ real(kind=8) :: phi
+ real(kind=8),dimension(0:2,0:5) :: R
+
+ locmove_output=.true.
+
+! call angles2tab(30.*deg2rad,70.*deg2rad,a_n,a_ang,a_tab)
+! call angles2tab(80.*deg2rad,130.*deg2rad,b_n,b_ang,b_tab)
+! call minmax_angles(0.D0,3.8D0,0.D0,3.8D0,b_n,b_ang,b_tab)
+! call construct_tab
+! call output_tabs
+
+! call construct_ranges(phi_n,phi_start,phi_end)
+! do i=0,phi_n-1
+! print *,phi_start(i)*rad2deg,phi_end(i)*rad2deg
+! enddo
+
+! call fix_no_moves(phi)
+! print *,'NO MOVES FOUND, BEST PHI IS',phi*rad2deg
+
+ R(0,0)=0.D0
+ R(1,0)=0.D0
+ R(2,0)=0.D0
+ R(0,1)=0.D0
+ R(1,1)=-cos(28.D0*deg2rad)
+ R(2,1)=-0.5D0-sin(28.D0*deg2rad)
+ R(0,2)=0.D0
+ R(1,2)=0.D0
+ R(2,2)=-0.5D0
+ R(0,3)=cos(30.D0*deg2rad)
+ R(1,3)=0.D0
+ R(2,3)=0.D0
+ R(0,4)=0.D0
+ R(1,4)=0.D0
+ R(2,4)=0.5D0
+ R(0,5)=0.D0
+ R(1,5)=cos(26.D0*deg2rad)
+ R(2,5)=0.5D0+sin(26.D0*deg2rad)
+ do i=1,5
+ do j=0,2
+ R(j,i)=vbl*R(j,i)
+ enddo
+ enddo
+! i=move_res(R(0,1),0.D0*deg2rad,180.D0*deg2rad)
+ imov=2
+ i=move_res(0.D0*deg2rad,180.D0*deg2rad,imov)
+ print *,'RETURNED ',i
+ print *,(R(i,3)/vbl,i=0,2)
+
+ return
+ end subroutine loc_test
+#endif
+!-----------------------------------------------------------------------------
+! matmult.f
+!-----------------------------------------------------------------------------
+ subroutine MATMULT(A1,A2,A3)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+!el local variables
+ integer :: i,j,k
+ real(kind=8) :: A3IJ
+
+ real(kind=8),DIMENSION(3,3) :: A1,A2,A3
+ real(kind=8),DIMENSION(3,3) :: AI3
+ DO 1 I=1,3
+ DO 2 J=1,3
+ A3IJ=0.0
+ DO 3 K=1,3
+ 3 A3IJ=A3IJ+A1(I,K)*A2(K,J)
+ AI3(I,J)=A3IJ
+ 2 CONTINUE
+ 1 CONTINUE
+ DO 4 I=1,3
+ DO 4 J=1,3
+ 4 A3(I,J)=AI3(I,J)
+ return
+ end subroutine MATMULT
+!-----------------------------------------------------------------------------
+! readpdb.F
+!-----------------------------------------------------------------------------
+ subroutine int_from_cart(lside,lprn)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+ use control_data,only:out1file
+#ifdef MPI
+ include "mpif.h"
+#endif
+! include 'COMMON.LOCAL'
+! include 'COMMON.VAR'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+! include 'COMMON.NAMES'
+! include 'COMMON.CONTROL'
+! include 'COMMON.SETUP'
+ character(len=3) :: seq,res
+! character*5 atom
+ character(len=80) :: card
+ real(kind=8),dimension(3,20) :: sccor
+ integer :: i,j,iti !el rescode,
+ logical :: lside,lprn
+ real(kind=8) :: di,cosfac,sinfac
+ integer :: nres2
+ nres2=2*nres
+
+ if(me.eq.king.or..not.out1file)then
+ if (lprn) then
+ write (iout,'(/a)') &
+ 'Internal coordinates calculated from crystal structure.'
+ if (lside) then
+ write (iout,'(8a)') ' Res ',' dvb',' Theta',&
+ ' Gamma',' Dsc_id',' Dsc',' Alpha',&
+ ' Beta '
+ else
+ write (iout,'(4a)') ' Res ',' dvb',' Theta',&
+ ' Gamma'
+ endif
+ endif
+ endif
+ do i=1,nres-1
+!in wham do i=1,nres
+ iti=itype(i)
+ if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
+ write (iout,'(a,i4)') 'Bad Cartesians for residue',i
+!test stop
+ endif
+!#ifndef WHAM_RUN
+ vbld(i+1)=dist(i,i+1)
+ vbld_inv(i+1)=1.0d0/vbld(i+1)
+!#endif
+ if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
+ if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
+ enddo
+!el -----
+!#ifdef WHAM_RUN
+! if (itype(1).eq.ntyp1) then
+! do j=1,3
+! c(j,1)=c(j,2)+(c(j,3)-c(j,4))
+! enddo
+! endif
+! if (itype(nres).eq.ntyp1) then
+! do j=1,3
+! c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3))
+! enddo
+! endif
+!#endif
+! if (unres_pdb) then
+! if (itype(1).eq.21) then
+! theta(3)=90.0d0*deg2rad
+! phi(4)=180.0d0*deg2rad
+! vbld(2)=3.8d0
+! vbld_inv(2)=1.0d0/vbld(2)
+! endif
+! if (itype(nres).eq.21) then
+! theta(nres)=90.0d0*deg2rad
+! phi(nres)=180.0d0*deg2rad
+! vbld(nres)=3.8d0
+! vbld_inv(nres)=1.0d0/vbld(2)
+! endif
+! endif
+ if (lside) then
+ do i=2,nres-1
+ do j=1,3
+ c(j,nres2+2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i) &
+ +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
+! in wham c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1)
+ enddo
+ iti=itype(i)
+ di=dist(i,nres+i)
+!#ifndef WHAM_RUN
+! 10/03/12 Adam: Correction for zero SC-SC bond length
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1 .and. di.eq.0.0d0) &
+ di=dsc(itype(i))
+ vbld(i+nres)=di
+ if (itype(i).ne.10) then
+ vbld_inv(i+nres)=1.0d0/di
+ else
+ vbld_inv(i+nres)=0.0d0
+ endif
+!#endif
+ if (iti.ne.10) then
+ alph(i)=alpha(nres+i,i,nres2+2)
+ omeg(i)=beta(nres+i,i,nres2+2,i+1)
+ endif
+ if(me.eq.king.or..not.out1file)then
+ if (lprn) &
+ write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),&
+ rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),&
+ rad2deg*alph(i),rad2deg*omeg(i)
+ endif
+ enddo
+ else if (lprn) then
+ do i=2,nres
+ iti=itype(i)
+ if(me.eq.king.or..not.out1file) &
+ write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),&
+ rad2deg*theta(i),rad2deg*phi(i)
+ enddo
+ endif
+ return
+ end subroutine int_from_cart
+!-----------------------------------------------------------------------------
+ subroutine sc_loc_geom(lprn)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+ use control_data,only:out1file
+#ifdef MPI
+ include "mpif.h"
+#endif
+! include 'COMMON.LOCAL'
+! include 'COMMON.VAR'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+! include 'COMMON.NAMES'
+! include 'COMMON.CONTROL'
+! include 'COMMON.SETUP'
+ real(kind=8),dimension(3) :: x_prime,y_prime,z_prime
+ logical :: lprn
+!el local variables
+ integer :: i,j,it,iti
+ real(kind=8) :: cosfac2,sinfac2,xx,yy,zz,cosfac,sinfac
+ do i=1,nres-1
+ do j=1,3
+ dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
+ enddo
+ enddo
+ do i=2,nres-1
+ if (itype(i).ne.10) then
+ do j=1,3
+ dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
+ enddo
+ else
+ do j=1,3
+ dc_norm(j,i+nres)=0.0d0
+ enddo
+ endif
+ enddo
+ do i=2,nres-1
+ costtab(i+1) =dcos(theta(i+1))
+ sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
+ cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
+ sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
+ cosfac2=0.5d0/(1.0d0+costtab(i+1))
+ cosfac=dsqrt(cosfac2)
+ sinfac2=0.5d0/(1.0d0-costtab(i+1))
+ sinfac=dsqrt(sinfac2)
+ it=itype(i)
+
+ if ((it.ne.10).and.(it.ne.ntyp1)) then
+!el if (it.ne.10) then
+!
+! Compute the axes of tghe local cartesian coordinates system; store in
+! x_prime, y_prime and z_prime
+!
+ do j=1,3
+ x_prime(j) = 0.00
+ y_prime(j) = 0.00
+ z_prime(j) = 0.00
+ enddo
+ do j = 1,3
+ x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
+ y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
+ enddo
+ call vecpr(x_prime,y_prime,z_prime)
+!
+! Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
+! to local coordinate system. Store in xx, yy, zz.
+!
+ xx=0.0d0
+ yy=0.0d0
+ zz=0.0d0
+ do j = 1,3
+ xx = xx + x_prime(j)*dc_norm(j,i+nres)
+ yy = yy + y_prime(j)*dc_norm(j,i+nres)
+ zz = zz + z_prime(j)*dc_norm(j,i+nres)
+ enddo
+
+ xxref(i)=xx
+ yyref(i)=yy
+ zzref(i)=zz
+ else
+ xxref(i)=0.0d0
+ yyref(i)=0.0d0
+ zzref(i)=0.0d0
+ endif
+ enddo
+ if (lprn) then
+ do i=2,nres
+ iti=itype(i)
+ if(me.eq.king.or..not.out1file) &
+ write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),&
+ yyref(i),zzref(i)
+ enddo
+ endif
+
+ return
+ end subroutine sc_loc_geom
+!-----------------------------------------------------------------------------
+ subroutine sccenter(ires,nscat,sccor)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CHAIN'
+ integer :: i,j,ires,nscat
+ real(kind=8),dimension(3,20) :: sccor
+ real(kind=8) :: sccmj
+ do j=1,3
+ sccmj=0.0D0
+ do i=1,nscat
+ sccmj=sccmj+sccor(j,i)
+ enddo
+ dc(j,ires)=sccmj/nscat
+ enddo
+ return
+ end subroutine sccenter
+#if !defined(WHAM_RUN) && !defined(CLUSTER)
+!-----------------------------------------------------------------------------
+ subroutine bond_regular
+ use calc_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.VAR'
+! include 'COMMON.LOCAL'
+! include 'COMMON.CALC'
+! include 'COMMON.INTERACT'
+! include 'COMMON.CHAIN'
+ do i=1,nres-1
+ vbld(i+1)=vbl
+ vbld_inv(i+1)=1.0d0/vbld(i+1)
+ vbld(i+1+nres)=dsc(itype(i+1))
+ vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
+! print *,vbld(i+1),vbld(i+1+nres)
+ enddo
+ return
+ end subroutine bond_regular
+#endif
+!-----------------------------------------------------------------------------
+! refsys.f
+!-----------------------------------------------------------------------------
+ subroutine refsys(i2,i3,i4,e1,e2,e3,fail)
+! This subroutine calculates unit vectors of a local reference system
+! defined by atoms (i2), (i3), and (i4). The x axis is the axis from
+! atom (i3) to atom (i2), and the xy plane is the plane defined by atoms
+! (i2), (i3), and (i4). z axis is directed according to the sign of the
+! vector product (i3)-(i2) and (i3)-(i4). Sets fail to .true. if atoms
+! (i2) and (i3) or (i3) and (i4) coincide or atoms (i2), (i3), and (i4)
+! form a linear fragment. Returns vectors e1, e2, and e3.
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+ logical :: fail
+ real(kind=8),dimension(3) :: e1,e2,e3
+ real(kind=8),dimension(3) :: u,z
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CHAIN'
+ real(kind=8) :: coinc=1.0D-13,align=1.0D-13
+!el local variables
+ integer :: i,i1,i2,i3,i4
+ real(kind=8) :: v1,v2,v3,s1,s2,zi,ui,anorm
+ fail=.false.
+ s1=0.0
+ s2=0.0
+ do 1 i=1,3
+ zi=c(i,i2)-c(i,i3)
+ ui=c(i,i4)-c(i,i3)
+ s1=s1+zi*zi
+ s2=s2+ui*ui
+ z(i)=zi
+ 1 u(i)=ui
+ s1=sqrt(s1)
+ s2=sqrt(s2)
+ if (s1.gt.coinc) goto 2
+ write (iout,1000) i2,i3,i1
+ fail=.true.
+! do 3 i=1,3
+! 3 c(i,i1)=0.0D0
+ return
+ 2 if (s2.gt.coinc) goto 4
+ write(iout,1000) i3,i4,i1
+ fail=.true.
+ do 5 i=1,3
+ 5 c(i,i1)=0.0D0
+ return
+ 4 s1=1.0/s1
+ s2=1.0/s2
+ v1=z(2)*u(3)-z(3)*u(2)
+ v2=z(3)*u(1)-z(1)*u(3)
+ v3=z(1)*u(2)-z(2)*u(1)
+ anorm=dsqrt(v1*v1+v2*v2+v3*v3)
+ if (anorm.gt.align) goto 6
+ write (iout,1010) i2,i3,i4,i1
+ fail=.true.
+! do 7 i=1,3
+! 7 c(i,i1)=0.0D0
+ return
+ 6 anorm=1.0D0/anorm
+ e3(1)=v1*anorm
+ e3(2)=v2*anorm
+ e3(3)=v3*anorm
+ e1(1)=z(1)*s1
+ e1(2)=z(2)*s1
+ e1(3)=z(3)*s1
+ e2(1)=e1(3)*e3(2)-e1(2)*e3(3)
+ e2(2)=e1(1)*e3(3)-e1(3)*e3(1)
+ e2(3)=e1(2)*e3(1)-e1(1)*e3(2)
+ 1000 format (/1x,' * * * Error - atoms',i4,' and',i4,' coincide.',&
+ 'coordinates of atom',i4,' are set to zero.')
+ 1010 format (/1x,' * * * Error - atoms',2(i4,2h, ),i4,' form a linear',&
+ ' fragment. coordinates of atom',i4,' are set to zero.')
+ return
+ end subroutine refsys
+!-----------------------------------------------------------------------------
+! int_to_cart.f
+!-----------------------------------------------------------------------------
+ subroutine int_to_cart
+!--------------------------------------------------------------
+! This subroutine converts the energy derivatives from internal
+! coordinates to cartesian coordinates
+!-------------------------------------------------------------
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.VAR'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.GEO'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.MD'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.SCCOR'
+! calculating dE/ddc1
+!el local variables
+ integer :: j,i
+ if (nres.lt.3) go to 18
+ do j=1,3
+ gcart(j,1)=gcart(j,1)+gloc(1,icg)*dphi(j,1,4) &
+ +gloc(nres-2,icg)*dtheta(j,1,3)
+ if(itype(2).ne.10) then
+ gcart(j,1)=gcart(j,1)+gloc(ialph(2,1),icg)*dalpha(j,1,2)+ &
+ gloc(ialph(2,1)+nside,icg)*domega(j,1,2)
+ endif
+ enddo
+! Calculating the remainder of dE/ddc2
+ do j=1,3
+ gcart(j,2)=gcart(j,2)+gloc(1,icg)*dphi(j,2,4)+ &
+ gloc(nres-2,icg)*dtheta(j,2,3)+gloc(nres-1,icg)*dtheta(j,1,4)
+ if(itype(2).ne.10) then
+ gcart(j,2)=gcart(j,2)+gloc(ialph(2,1),icg)*dalpha(j,2,2)+ &
+ gloc(ialph(2,1)+nside,icg)*domega(j,2,2)
+ endif
+ if(itype(3).ne.10) then
+ gcart(j,2)=gcart(j,2)+gloc(ialph(3,1),icg)*dalpha(j,1,3)+ &
+ gloc(ialph(3,1)+nside,icg)*domega(j,1,3)
+ endif
+ if(nres.gt.4) then
+ gcart(j,2)=gcart(j,2)+gloc(2,icg)*dphi(j,1,5)
+ endif
+ enddo
+! If there are only five residues
+ if(nres.eq.5) then
+ do j=1,3
+ gcart(j,3)=gcart(j,3)+gloc(1,icg)*dphi(j,3,4)+gloc(2,icg)* &
+ dphi(j,2,5)+gloc(nres-1,icg)*dtheta(j,2,4)+gloc(nres,icg)* &
+ dtheta(j,1,5)
+ if(itype(3).ne.10) then
+ gcart(j,3)=gcart(j,3)+gloc(ialph(3,1),icg)* &
+ dalpha(j,2,3)+gloc(ialph(3,1)+nside,icg)*domega(j,2,3)
+ endif
+ if(itype(4).ne.10) then
+ gcart(j,3)=gcart(j,3)+gloc(ialph(4,1),icg)* &
+ dalpha(j,1,4)+gloc(ialph(4,1)+nside,icg)*domega(j,1,4)
+ endif
+ enddo
+ endif
+! If there are more than five residues
+ if(nres.gt.5) then
+ do i=3,nres-3
+ do j=1,3
+ gcart(j,i)=gcart(j,i)+gloc(i-2,icg)*dphi(j,3,i+1) &
+ +gloc(i-1,icg)*dphi(j,2,i+2)+ &
+ gloc(i,icg)*dphi(j,1,i+3)+gloc(nres+i-4,icg)*dtheta(j,2,i+1)+ &
+ gloc(nres+i-3,icg)*dtheta(j,1,i+2)
+ if(itype(i).ne.10) then
+ gcart(j,i)=gcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,2,i)+ &
+ gloc(ialph(i,1)+nside,icg)*domega(j,2,i)
+ endif
+ if(itype(i+1).ne.10) then
+ gcart(j,i)=gcart(j,i)+gloc(ialph(i+1,1),icg)*dalpha(j,1,i+1) &
+ +gloc(ialph(i+1,1)+nside,icg)*domega(j,1,i+1)
+ endif
+ enddo
+ enddo
+ endif
+! Setting dE/ddnres-2
+ if(nres.gt.5) then
+ do j=1,3
+ gcart(j,nres-2)=gcart(j,nres-2)+gloc(nres-4,icg)* &
+ dphi(j,3,nres-1)+gloc(nres-3,icg)*dphi(j,2,nres) &
+ +gloc(2*nres-6,icg)* &
+ dtheta(j,2,nres-1)+gloc(2*nres-5,icg)*dtheta(j,1,nres)
+ if(itype(nres-2).ne.10) then
+ gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-2,1),icg)* &
+ dalpha(j,2,nres-2)+gloc(ialph(nres-2,1)+nside,icg)* &
+ domega(j,2,nres-2)
+ endif
+ if(itype(nres-1).ne.10) then
+ gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-1,1),icg)* &
+ dalpha(j,1,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* &
+ domega(j,1,nres-1)
+ endif
+ enddo
+ endif
+! Settind dE/ddnres-1
+ do j=1,3
+ gcart(j,nres-1)=gcart(j,nres-1)+gloc(nres-3,icg)*dphi(j,3,nres)+ &
+ gloc(2*nres-5,icg)*dtheta(j,2,nres)
+ if(itype(nres-1).ne.10) then
+ gcart(j,nres-1)=gcart(j,nres-1)+gloc(ialph(nres-1,1),icg)* &
+ dalpha(j,2,nres-1)+gloc(ialph(nres-1,1)+nside,icg)* &
+ domega(j,2,nres-1)
+ endif
+ enddo
+! The side-chain vector derivatives
+ do i=2,nres-1
+ if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ do j=1,3
+ gxcart(j,i)=gxcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,3,i) &
+ +gloc(ialph(i,1)+nside,icg)*domega(j,3,i)
+ enddo
+ endif
+ enddo
+!----------------------------------------------------------------------
+! INTERTYP=1 SC...Ca...Ca...Ca
+! INTERTYP=2 Ca...Ca...Ca...SC
+! INTERTYP=3 SC...Ca...Ca...SC
+! calculating dE/ddc1
+ 18 continue
+! do i=1,nres
+! gloc(i,icg)=0.0D0
+! write (iout,*) "poczotkoawy",i,gloc_sc(1,i,icg)
+! enddo
+ if (nres.lt.2) return
+ if ((nres.lt.3).and.(itype(1).eq.10)) return
+ if ((itype(1).ne.10).and.(itype(1).ne.ntyp1)) then
+ do j=1,3
+!c Derviative was calculated for oposite vector of side chain therefore
+! there is "-" sign before gloc_sc
+ gxcart(j,1)=gxcart(j,1)-gloc_sc(1,0,icg)* &
+ dtauangle(j,1,1,3)
+ gcart(j,1)=gcart(j,1)+gloc_sc(1,0,icg)* &
+ dtauangle(j,1,2,3)
+ if ((itype(2).ne.10).and.(itype(2).ne.ntyp1)) then
+ gxcart(j,1)= gxcart(j,1) &
+ -gloc_sc(3,0,icg)*dtauangle(j,3,1,3)
+ gcart(j,1)=gcart(j,1)+gloc_sc(3,0,icg)* &
+ dtauangle(j,3,2,3)
+ endif
+ enddo
+ endif
+ if ((nres.ge.3).and.(itype(3).ne.10).and.(itype(3).ne.ntyp1)) &
+ then
+ do j=1,3
+ gcart(j,1)=gcart(j,1)+gloc_sc(2,1,icg)*dtauangle(j,2,1,4)
+ enddo
+ endif
+! As potetnial DO NOT depend on omicron anlge their derivative is
+! ommited
+! & +gloc_sc(intertyp,nres-2,icg)*dtheta(j,1,3)
+
+! Calculating the remainder of dE/ddc2
+ do j=1,3
+ if((itype(2).ne.10).and.(itype(2).ne.ntyp1)) then
+ if (itype(1).ne.10) gxcart(j,2)=gxcart(j,2)+ &
+ gloc_sc(3,0,icg)*dtauangle(j,3,3,3)
+ if ((itype(3).ne.10).and.(nres.ge.3).and.(itype(3).ne.ntyp1)) &
+ then
+ gxcart(j,2)=gxcart(j,2)-gloc_sc(3,1,icg)*dtauangle(j,3,1,4)
+!c the - above is due to different vector direction
+ gcart(j,2)=gcart(j,2)+gloc_sc(3,1,icg)*dtauangle(j,3,2,4)
+ endif
+ if (nres.gt.3) then
+ gxcart(j,2)=gxcart(j,2)-gloc_sc(1,1,icg)*dtauangle(j,1,1,4)
+!c the - above is due to different vector direction
+ gcart(j,2)=gcart(j,2)+gloc_sc(1,1,icg)*dtauangle(j,1,2,4)
+! write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,2,4),"gcart"
+! write(iout,*) gloc_sc(1,1,icg),dtauangle(j,1,1,4),"gx"
+ endif
+ endif
+ if ((itype(1).ne.10).and.(itype(1).ne.ntyp1)) then
+ gcart(j,2)=gcart(j,2)+gloc_sc(1,0,icg)*dtauangle(j,1,3,3)
+! write(iout,*) gloc_sc(1,0,icg),dtauangle(j,1,3,3)
+ endif
+ if ((itype(3).ne.10).and.(nres.ge.3)) then
+ gcart(j,2)=gcart(j,2)+gloc_sc(2,1,icg)*dtauangle(j,2,2,4)
+! write(iout,*) gloc_sc(2,1,icg),dtauangle(j,2,2,4)
+ endif
+ if ((itype(4).ne.10).and.(nres.ge.4)) then
+ gcart(j,2)=gcart(j,2)+gloc_sc(2,2,icg)*dtauangle(j,2,1,5)
+! write(iout,*) gloc_sc(2,2,icg),dtauangle(j,2,1,5)
+ endif
+
+! write(iout,*) gcart(j,2),itype(2),itype(1),itype(3), "gcart2"
+ enddo
+! If there are more than five residues
+ if(nres.ge.5) then
+ do i=3,nres-2
+ do j=1,3
+! write(iout,*) "before", gcart(j,i)
+ if ((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
+ gxcart(j,i)=gxcart(j,i)+gloc_sc(2,i-2,icg) &
+ *dtauangle(j,2,3,i+1) &
+ -gloc_sc(1,i-1,icg)*dtauangle(j,1,1,i+2)
+ gcart(j,i)=gcart(j,i)+gloc_sc(1,i-1,icg) &
+ *dtauangle(j,1,2,i+2)
+! write(iout,*) "new",j,i,
+! & gcart(j,i),gloc_sc(1,i-1,icg),dtauangle(j,1,2,i+2)
+ if (itype(i-1).ne.10) then
+ gxcart(j,i)=gxcart(j,i)+gloc_sc(3,i-2,icg) &
+ *dtauangle(j,3,3,i+1)
+ endif
+ if (itype(i+1).ne.10) then
+ gxcart(j,i)=gxcart(j,i)-gloc_sc(3,i-1,icg) &
+ *dtauangle(j,3,1,i+2)
+ gcart(j,i)=gcart(j,i)+gloc_sc(3,i-1,icg) &
+ *dtauangle(j,3,2,i+2)
+ endif
+ endif
+ if (itype(i-1).ne.10) then
+ gcart(j,i)=gcart(j,i)+gloc_sc(1,i-2,icg)* &
+ dtauangle(j,1,3,i+1)
+ endif
+ if (itype(i+1).ne.10) then
+ gcart(j,i)=gcart(j,i)+gloc_sc(2,i-1,icg)* &
+ dtauangle(j,2,2,i+2)
+! write(iout,*) "numer",i,gloc_sc(2,i-1,icg),
+! & dtauangle(j,2,2,i+2)
+ endif
+ if (itype(i+2).ne.10) then
+ gcart(j,i)=gcart(j,i)+gloc_sc(2,i,icg)* &
+ dtauangle(j,2,1,i+3)
+ endif
+ enddo
+ enddo
+ endif
+! Setting dE/ddnres-1
+ if(nres.ge.4) then
+ do j=1,3
+ if ((itype(nres-1).ne.10).and.(itype(nres-1).ne.ntyp1)) then
+ gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(2,nres-3,icg) &
+ *dtauangle(j,2,3,nres)
+! write (iout,*) "gxcart(nres-1)", gloc_sc(2,nres-3,icg),
+! & dtauangle(j,2,3,nres), gxcart(j,nres-1)
+ if (itype(nres-2).ne.10) then
+ gxcart(j,nres-1)=gxcart(j,nres-1)+gloc_sc(3,nres-3,icg) &
+ *dtauangle(j,3,3,nres)
+ endif
+ if ((itype(nres).ne.10).and.(itype(nres).ne.ntyp1)) then
+ gxcart(j,nres-1)=gxcart(j,nres-1)-gloc_sc(3,nres-2,icg) &
+ *dtauangle(j,3,1,nres+1)
+ gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(3,nres-2,icg) &
+ *dtauangle(j,3,2,nres+1)
+ endif
+ endif
+ if ((itype(nres-2).ne.10).and.(itype(nres-2).ne.ntyp1)) then
+ gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(1,nres-3,icg)* &
+ dtauangle(j,1,3,nres)
+ endif
+ if ((itype(nres).ne.10).and.(itype(nres).ne.ntyp1)) then
+ gcart(j,nres-1)=gcart(j,nres-1)+gloc_sc(2,nres-2,icg)* &
+ dtauangle(j,2,2,nres+1)
+! write (iout,*) "gcart(nres-1)", gloc_sc(2,nres-2,icg),
+! & dtauangle(j,2,2,nres+1), itype(nres-1),itype(nres)
+ endif
+ enddo
+ endif
+! Settind dE/ddnres
+ if ((nres.ge.3).and.(itype(nres).ne.10).and. &
+ (itype(nres).ne.ntyp1))then
+ do j=1,3
+ gxcart(j,nres)=gxcart(j,nres)+gloc_sc(3,nres-2,icg) &
+ *dtauangle(j,3,3,nres+1)+gloc_sc(2,nres-2,icg) &
+ *dtauangle(j,2,3,nres+1)
+ enddo
+ endif
+! The side-chain vector derivatives
+ return
+ end subroutine int_to_cart
+#if !defined(WHAM_RUN) && !defined(CLUSTER)
+!-----------------------------------------------------------------------------
+! readrtns_CSA.F
+!-----------------------------------------------------------------------------
+ subroutine gen_dist_constr
+! Generate CA distance constraints.
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.INTERACT'
+! include 'COMMON.LOCAL'
+! include 'COMMON.NAMES'
+! include 'COMMON.CHAIN'
+! include 'COMMON.FFIELD'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.HEADER'
+! include 'COMMON.CONTROL'
+! include 'COMMON.DBASE'
+! include 'COMMON.THREAD'
+! include 'COMMON.TIME1'
+! integer :: itype_pdb !(maxres)
+! common /pizda/ itype_pdb(nres)
+ character(len=2) :: iden
+!el local variables
+ integer :: i,j
+!d print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
+!d write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
+!d & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
+!d & ' nsup',nsup
+ do i=nstart_sup,nstart_sup+nsup-1
+!d write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
+!d & ' seq_pdb', restyp(itype_pdb(i))
+ do j=i+2,nstart_sup+nsup-1
+ nhpb=nhpb+1
+ ihpb(nhpb)=i+nstart_seq-nstart_sup
+ jhpb(nhpb)=j+nstart_seq-nstart_sup
+ forcon(nhpb)=weidis
+ dhpb(nhpb)=dist(i,j)
+ enddo
+ enddo
+!d write (iout,'(a)') 'Distance constraints:'
+!d do i=nss+1,nhpb
+!d ii=ihpb(i)
+!d jj=jhpb(i)
+!d iden='CA'
+!d if (ii.gt.nres) then
+!d iden='SC'
+!d ii=ii-nres
+!d jj=jj-nres
+!d endif
+!d write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
+!d & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
+!d & dhpb(i),forcon(i)
+!d enddo
+! deallocate(itype_pdb)
+
+ return
+ end subroutine gen_dist_constr
+#endif
+!-----------------------------------------------------------------------------
+! cartprint.f
+!-----------------------------------------------------------------------------
+ subroutine cartprint
+
+ use geometry_data, only: c
+ use energy_data, only: itype
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.NAMES'
+! include 'COMMON.IOUNITS'
+ integer :: i
+
+ write (iout,100)
+ do i=1,nres
+ write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),&
+ c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
+ enddo
+ 100 format (//' alpha-carbon coordinates ',&
+ ' centroid coordinates'/ &
+ ' ', 6X,'X',11X,'Y',11X,'Z',&
+ 10X,'X',11X,'Y',11X,'Z')
+ 110 format (a,'(',i3,')',6f12.5)
+ return
+ end subroutine cartprint
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+ subroutine alloc_geo_arrays
+!EL Allocation of tables used by module energy
+
+ integer :: i,j,nres2
+ nres2=2*nres
+! commom.bounds
+! common /bounds/
+ allocate(phibound(2,nres+2)) !(2,maxres)
+!----------------------
+! commom.chain
+! common /chain/ in molread
+! real(kind=8),dimension(:,:),allocatable :: c !(3,maxres2+2)
+! real(kind=8),dimension(:,:),allocatable :: dc
+ allocate(dc_old(3,0:nres2))
+! if(.not.allocated(dc_norm2)) allocate(dc_norm2(3,0:nres2+2)) !(3,0:maxres2)
+ if(.not.allocated(dc_norm2)) then
+ allocate(dc_norm2(3,0:nres2+2)) !(3,0:maxres2)
+ dc_norm2(:,:)=0.d0
+ endif
+!
+!el if(.not.allocated(dc_norm))
+!elwrite(iout,*) "jestem w alloc geo 1"
+ if(.not.allocated(dc_norm)) then
+ allocate(dc_norm(3,0:nres2+2)) !(3,0:maxres2)
+ dc_norm(:,:)=0.d0
+ endif
+!elwrite(iout,*) "jestem w alloc geo 1"
+ allocate(xloc(3,nres),xrot(3,nres))
+!elwrite(iout,*) "jestem w alloc geo 1"
+ xloc(:,:)=0.0D0
+!elwrite(iout,*) "jestem w alloc geo 1"
+ allocate(dc_work(6*nres)) !(MAXRES6) maxres6=6*maxres
+! common /rotmat/
+ allocate(t(3,3,nres),r(3,3,nres))
+ allocate(prod(3,3,nres),rt(3,3,nres)) !(3,3,maxres)
+! common /refstruct/
+ if(.not.allocated(cref)) allocate(cref(3,nres2+2,maxperm)) !(3,maxres2+2,maxperm)
+!elwrite(iout,*) "jestem w alloc geo 2"
+ allocate(crefjlee(3,nres2+2)) !(3,maxres2+2)
+ if(.not.allocated(chain_rep)) allocate(chain_rep(3,nres2+2,maxsym)) !(3,maxres2+2,maxsym)
+ if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
+! common /from_zscore/ in module.compare
+!----------------------
+! common.local
+! Inverses of the actual virtual bond lengths
+! common /invlen/ in io_conf: molread or readpdb
+! real(kind=8),dimension(:),allocatable :: vbld_inv !(maxres2)
+!----------------------
+! common.var
+! Store the geometric variables in the following COMMON block.
+! common /var/ in readpdb or ...
+ if(.not.allocated(theta)) allocate(theta(nres+2))
+ if(.not.allocated(phi)) allocate(phi(nres+2))
+ if(.not.allocated(alph)) allocate(alph(nres+2))
+ if(.not.allocated(omeg)) allocate(omeg(nres+2))
+ if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
+ if(.not.allocated(phiref)) allocate(phiref(nres+2))
+ if(.not.allocated(costtab)) allocate(costtab(nres))
+ if(.not.allocated(sinttab)) allocate(sinttab(nres))
+ if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
+ if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
+! real(kind=8),dimension(:),allocatable :: vbld !(2*maxres) in io_conf: molread or readpdb
+ allocate(omicron(2,nres+2)) !(2,maxres)
+ allocate(tauangle(3,nres+2)) !(3,maxres)
+!elwrite(iout,*) "jestem w alloc geo 3"
+ if(.not.allocated(xxtab)) allocate(xxtab(nres))
+ if(.not.allocated(yytab)) allocate(yytab(nres))
+ if(.not.allocated(zztab)) allocate(zztab(nres)) !(maxres)
+ if(.not.allocated(xxref)) allocate(xxref(nres))
+ if(.not.allocated(yyref)) allocate(yyref(nres))
+ if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
+ allocate(ialph(nres,2)) !(maxres,2)
+ ialph(:,1)=0
+ ialph(:,2)=0
+ allocate(ivar(4*nres2)) !(4*maxres2)
+
+#if defined(WHAM_RUN) || defined(CLUSTER)
+ allocate(vbld(2*nres))
+ vbld(:)=0.d0
+ allocate(vbld_inv(2*nres))
+ vbld_inv(:)=0.d0
+#endif
+
+ return
+ end subroutine alloc_geo_arrays
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+ end module geometry
projects / unres4.git / blobdiff