diagnoza dlaczego WHAM nie dziala
[unres.git] / source / wham / src-M / energy_p_new.F
index a4f5fb4..810275a 100644 (file)
@@ -44,11 +44,13 @@ C Gay-Berne potential (shifted LJ, angular dependence).
       goto 106
 C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence).
   105 call egbv(evdw,evdw_t)
+C      write(iout,*) 'po elektostatyce'
 C
 C Calculate electrostatic (H-bonding) energy of the main chain.
 C
-  106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
-C
+  106  call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
+C            write(iout,*) 'po eelec'
+
 C Calculate excluded-volume interaction energy between peptide groups
 C and side chains.
 C
@@ -56,8 +58,9 @@ C
 c
 c Calculate the bond-stretching energy
 c
+
       call ebond(estr)
-c      write (iout,*) "estr",estr
+C       write (iout,*) "estr",estr
 C 
 C Calculate the disulfide-bridge and other energy and the contributions
 C from other distance constraints.
@@ -68,12 +71,12 @@ C
 C Calculate the virtual-bond-angle energy.
 C
       call ebend(ebe)
-cd    print *,'Bend energy finished.'
+C      print *,'Bend energy finished.'
 C
 C Calculate the SC local energy.
 C
       call esc(escloc)
-cd    print *,'SCLOC energy finished.'
+C       print *,'SCLOC energy finished.'
 C
 C Calculate the virtual-bond torsional energy.
 C
@@ -1898,13 +1901,14 @@ cd      enddo
       do i=1,nres
         num_cont_hb(i)=0
       enddo
-cd      print '(a)','Enter EELEC'
-cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
+C      print '(a)','Enter EELEC'
+C      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
       do i=1,nres
         gel_loc_loc(i)=0.0d0
         gcorr_loc(i)=0.0d0
       enddo
       do i=iatel_s,iatel_e
+          if (i.eq.1) cycle
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
      &  .or. itype(i+2).eq.ntyp1
      &  .or. itype(i-1).eq.ntyp1
@@ -1926,8 +1930,9 @@ cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
           zmedi=mod(zmedi,boxzsize)
           if (zmedi.lt.0) zmedi=zmedi+boxzsize
         num_conti=0
-c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
+C        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         do j=ielstart(i),ielend(i)
+          if (j.eq.1) cycle
           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
      & .or.itype(j+2).eq.ntyp1
      & .or.itype(j-1).eq.ntyp1
@@ -2027,7 +2032,7 @@ c             write (iout,'(a6,2i5,0pf7.3,2i5,2e11.3)')
 c     &'evdw1',i,j,evdwij
 c     &,iteli,itelj,aaa,evdw1
 
-c              write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
+              write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
 c          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
 c     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
 c     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
@@ -3272,7 +3277,8 @@ c          write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff
             gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
           enddo
         endif
-
+        write (iout,'(a7,i5,4f7.3)')
+     &     "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
       enddo
       estr=0.5d0*AKP*estr+estr1
 c
@@ -3357,6 +3363,7 @@ c     write (*,'(a,i2)') 'EBEND ICG=',icg
 c      write (iout,*) ithet_start,ithet_end
       do i=ithet_start,ithet_end
 C        if (itype(i-1).eq.ntyp1) cycle
+        if (i.le.2) cycle
         if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
      &  .or.itype(i).eq.ntyp1) cycle
 C Zero the energy function and its derivative at 0 or pi.
@@ -3374,6 +3381,10 @@ C Zero the energy function and its derivative at 0 or pi.
           ichir21=isign(1,itype(i))
           ichir22=isign(1,itype(i))
          endif
+         if (i.eq.3) then
+          y(1)=0.0D0
+          y(2)=0.0D0
+          else
 
         if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
 #ifdef OSF
@@ -3390,6 +3401,7 @@ C Zero the energy function and its derivative at 0 or pi.
           y(1)=0.0D0
           y(2)=0.0D0
         endif
+        endif
         if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
 #ifdef OSF
           phii1=phi(i+1)
@@ -3457,6 +3469,8 @@ C Derivatives of the "mean" values in gamma1 and gamma2.
      &        E_theta,E_tc)
         endif
         etheta=etheta+ethetai
+c         write (iout,'(a6,i5,0pf7.3,f7.3,i5)')
+c     &      'ebend',i,ethetai,theta(i),itype(i)
 c        write (iout,'(2i3,3f8.3,f10.5)') i,it,rad2deg*theta(i),
 c     &    rad2deg*phii,rad2deg*phii1,ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
@@ -3605,7 +3619,9 @@ C
       etheta=0.0D0
 c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
       do i=ithet_start,ithet_end
+C         if (i.eq.2) cycle
 C        if (itype(i-1).eq.ntyp1) cycle
+        if (i.le.2) cycle
         if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
      &  .or.itype(i).eq.ntyp1) cycle
         if (iabs(itype(i+1)).eq.20) iblock=2
@@ -3619,6 +3635,14 @@ C        if (itype(i-1).eq.ntyp1) cycle
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
         enddo
+        if (i.eq.3) then 
+          phii=0.0d0
+          ityp1=nthetyp+1
+          do k=1,nsingle
+            cosph1(k)=0.0d0
+            sinph1(k)=0.0d0
+          enddo
+        else
         if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
@@ -3639,6 +3663,7 @@ C        if (itype(i-1).eq.ntyp1) cycle
             sinph1(k)=0.0d0
           enddo 
         endif
+        endif
         if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
 #ifdef OSF
           phii1=phi(i+1)
@@ -3799,14 +3824,14 @@ C ALPHA and OMEGA.
       common /sccalc/ time11,time12,time112,theti,it,nlobit
       delta=0.02d0*pi
       escloc=0.0D0
-c     write (iout,'(a)') 'ESC'
+C      write (iout,*) 'ESC'
       do i=loc_start,loc_end
         it=itype(i)
         if (it.eq.ntyp1) cycle
         if (it.eq.10) goto 1
         nlobit=nlob(iabs(it))
 c       print *,'i=',i,' it=',it,' nlobit=',nlobit
-c       write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
+C        write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
         theti=theta(i+1)-pipol
         x(1)=dtan(theti)
         x(2)=alph(i)
@@ -3842,8 +3867,8 @@ c        write (iout,*) "i",i," x",x(1),x(2),x(3)
             dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
           enddo
           dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
-c         write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
-c    &             esclocbi,ss,ssd
+          write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
+     &             esclocbi,ss,ssd
           escloci=ss*escloci+(1.0d0-ss)*esclocbi
 c         escloci=esclocbi
 c         write (iout,*) escloci
@@ -3877,15 +3902,17 @@ c         write (iout,*) escloci
           enddo
           dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
 c         write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
-c    &             esclocbi,ss,ssd
+c     &             esclocbi,ss,ssd
           escloci=ss*escloci+(1.0d0-ss)*esclocbi
-c         write (iout,*) escloci
+C         write (iout,*) 'i=',i, escloci
         else
           call enesc(x,escloci,dersc,ddummy,.false.)
         endif
 
         escloc=escloc+escloci
-c        write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
+C        write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
+            write (iout,'(a6,i5,0pf7.3)')
+     &     'escloc',i,escloci
 
         gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
      &   wscloc*dersc(1)
@@ -4575,6 +4602,7 @@ C Set lprn=.true. for debugging
 c      lprn=.true.
       etors=0.0D0
       do i=iphi_start,iphi_end
+        if (i.le.2) cycle
         if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
      &      .or. itype(i).eq.ntyp1 .or. itype(i-3).eq.ntyp1) cycle
 C        if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
@@ -4678,6 +4706,7 @@ C Set lprn=.true. for debugging
 c     lprn=.true.
       etors_d=0.0D0
       do i=iphi_start,iphi_end-1
+        if (i.le.3) cycle
 C        if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1
 C     &      .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
          if ((itype(i-2).eq.ntyp1).or.itype(i-3).eq.ntyp1.or.