X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Funres%2Fgeometry.F90;fp=source%2Funres%2Fgeometry.F90;h=91cfcd41f4b9e313b1f09e353288d4095a04338a;hb=2611e37f82b576a1366c2d78fce87c1a55852037;hp=0000000000000000000000000000000000000000;hpb=a1e9d5a6b704c90ebd10f86a655780212a880dce;p=unres4.git diff --git a/source/unres/geometry.F90 b/source/unres/geometry.F90 new file mode 100644 index 0000000..91cfcd4 --- /dev/null +++ b/source/unres/geometry.F90 @@ -0,0 +1,3597 @@ + 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