ifdef poporawa
[unres4.git] / source / unres / geometry.f90
index 780684e..9737d1a 100644 (file)
@@ -9,6 +9,9 @@
       use energy_data
       implicit none
 !-----------------------------------------------------------------------------
+! commom.bounds
+!      common /bounds/
+!-----------------------------------------------------------------------------
 ! commom.chain
 !      common /chain/
 !      common /rotmat/
       nres2=2*nres
 ! Set lprn=.true. for debugging
       lprn = .false.
+      print *,"I ENTER CHAINBUILD"
 !
 ! 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.
 !
       write (iout,'(/a)') 'Recalculated internal coordinates'
       do i=2,nres-1
        do j=1,3
-         c(j,nres2)=0.5D0*(c(j,i-1)+c(j,i+1))  !maxres2=2*maxres
+         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,i+1)
+        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),be1
+        alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,nres2+2),be1
       enddo   
  1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
 
 
 #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
         dnorm1=dist(i-1,i)
         dnorm2=dist(i,i+1) 
        do j=1,3
-         c(j,nres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))/dnorm1 &
+         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
          tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
         endif
         endif
-        omeg(i)=beta(nres+i,i,nres2,i+1)
-        alph(i)=alpha(nres+i,i,nres2)
+        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)
 #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
 #endif
       return
       end subroutine int_from_cart1
-#ifndef WHAM_RUN
+#if !defined(WHAM_RUN) && !defined(CLUSTER)
 !-----------------------------------------------------------------------------
 ! check_sc_distr.f
 !-----------------------------------------------------------------------------
           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
       thetnorm=xx 
       return
       end function thetnorm
-#ifndef WHAM_RUN
+#if !defined(WHAM_RUN) && !defined(CLUSTER)
 !-----------------------------------------------------------------------------
       subroutine var_to_geom_restr(n,xx)
 !
 ! SCs.
       iteli=itel(i)
       do j=1,3
-       c(j,nres2+1)=0.5D0*(c(j,i)+c(j,i+1))
+!      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+1).lt.4.0D0*redfac) then
+       if (dist(nres+j,nres2+3).lt.4.0D0*redfac) then
          overlap=.true.
          return
         endif
 ! groups.
       do j=1,nnt-2
        do k=1,3
-         c(k,nres2+1)=0.5D0*(c(k,j)+c(k,j+1))
+         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+1).lt.4.0D0*redfac) then
+       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+2)=0.5D0*(c(j,i)+c(j,i+1))
+       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+2)=0.5D0*(c(k,j)+c(k,j+1))
+         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+1,nres2+2).lt.rpp(iteli,itelj)*redfac) then
+        if (dist(nres2+3,nres2+4).lt.rpp(iteli,itelj)*redfac) then
           overlap=.true.
           return
         endif
       endif
       if (box(1,2).lt.-MaxBoxLen .or. box(2,2).gt.MaxBoxLen) then
 #ifdef MPI
-        write (iout,'(a,i4,a)') 'CG Processor:',me,': bad sampling box.'
+        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.'
       dist=dsqrt(x12*x12+y12*y12+z12*z12)
       return
       end function dist
-#ifndef WHAM_RUN
+#if !defined(WHAM_RUN) && !defined(CLUSTER)
 !-----------------------------------------------------------------------------
 ! local_move.f
 !-----------------------------------------------------------------------------
        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
+        if ((dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0).and.&
+       (iti.ne.ntyp1  .and. itype(i+1).ne.ntyp1)) 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
       if (lside) then
         do i=2,nres-1
           do j=1,3
-            c(j,nres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i) &
+            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))
           else
             vbld_inv(i+nres)=0.0d0
           endif
+!#endif
           if (iti.ne.10) then
-            alph(i)=alpha(nres+i,i,nres2)
-            omeg(i)=beta(nres+i,i,nres2,i+1)
+            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) &
         sinfac2=0.5d0/(1.0d0-costtab(i+1))
         sinfac=dsqrt(sinfac2)
         it=itype(i)
-        if (it.ne.10) then
+
+        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 
             yyref(i),zzref(i)
         enddo
       endif
       return
       end subroutine sc_loc_geom
 !-----------------------------------------------------------------------------
       enddo
       return
       end subroutine sccenter
-#ifndef WHAM_RUN
+#if !defined(WHAM_RUN) && !defined(CLUSTER)
 !-----------------------------------------------------------------------------
       subroutine bond_regular
       use calc_data
 !   The side-chain vector derivatives
       return
       end subroutine int_to_cart
-#ifndef WHAM_RUN
+#if !defined(WHAM_RUN) && !defined(CLUSTER)
 !-----------------------------------------------------------------------------
 ! readrtns_CSA.F
 !-----------------------------------------------------------------------------
 !d     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
 !d     &  dhpb(i),forcon(i)
 !d      enddo
-      deallocate(itype_pdb)
+!      deallocate(itype_pdb)
 
       return
       end subroutine gen_dist_constr
 !      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)) !(3,0:maxres2)      
+!      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)) allocate(dc_norm(3,0:nres2)) !(3,0:maxres2)
+      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"
-      do i=1,nres
-       do j=1,3
-         xloc(j,i)=0.0D0
-        enddo
-      enddo
+      xloc(:,:)=0.0D0
 !elwrite(iout,*) "jestem w alloc geo 1"
       allocate(dc_work(6*nres)) !(MAXRES6) maxres6=6*maxres
 !      common /rotmat/
       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
 !-----------------------------------------------------------------------------