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