energi calculation working (no weights!)
[unres4.git] / source / unres / geometry.f90
index 7b6c961..e0b37d3 100644 (file)
       end subroutine sc_angular
 !-----------------------------------------------------------------------------
 ! initialize_p.F
+      subroutine sc_angular_nucl
+! 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'
+      use comm_locel
+      use calc_data_nucl
+      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
+      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_nucl
+
 !-----------------------------------------------------------------------------
       subroutine int_bounds(total_ints,lower_bound,upper_bound)
 !      implicit real*8 (a-h,o-z)
        endif
       endif
       do i=1,nres-1
+       if (molnum(i).ne.1) cycle
 !in wham      do i=1,nres
         iti=itype(i,1)
         if ((dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0).and.&
           di=dist(i,nres+i)
 !#ifndef WHAM_RUN
 ! 10/03/12 Adam: Correction for zero SC-SC bond length
+          
           if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1 .and. di.eq.0.0d0) &
-           di=dsc(itype(i,1))
+           di=dsc(itype(i,molnum(i)))
           vbld(i+nres)=di
           if (itype(i,1).ne.10) then
             vbld_inv(i+nres)=1.0d0/di
             alph(i)=alpha(nres+i,i,nres2+2)
             omeg(i)=beta(nres+i,i,nres2+2,i+1)
           endif
+          if (iti.ne.0) then
           if(me.eq.king.or..not.out1file)then
            if (lprn) &
            write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,vbld(i),&
            rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),&
            rad2deg*alph(i),rad2deg*omeg(i)
           endif
+          else
+          if(me.eq.king.or..not.out1file)then
+           if (lprn) &
+           write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,vbld(i),&
+           rad2deg*theta(i),rad2deg*phi(i),dsc(iti+1),vbld(nres+i),&
+           rad2deg*alph(i),rad2deg*omeg(i)
+          endif
+          endif
         enddo
       else if (lprn) then
         do i=2,nres
       integer :: i,j,ires,nscat
       real(kind=8),dimension(3,20) :: sccor
       real(kind=8) :: sccmj
-        print *,"I am in sccenter",ires,nscat
+!        print *,"I am in sccenter",ires,nscat
       do j=1,3
         sccmj=0.0D0
         do i=1,nscat
-          sccmj=sccmj+sccor(j,i) 
+          sccmj=sccmj+sccor(j,i)
+!C          print *,"insccent", ires,sccor(j,i) 
         enddo
         dc(j,ires)=sccmj/nscat
       enddo