criticall bug fix in energy

diff --git a/source/unres/energy.F90 b/source/unres/energy.F90
index 4e043fe..3870bb0 100644 (file)
--- a/source/unres/energy.F90
+++ b/source/unres/energy.F90
@@ -294,7 +294,7 @@
 !          print *,"Processor",myrank," BROADCAST iorder"
 ! FG master sets up the WEIGHTS_ array which will be broadcast to the 
 ! FG slaves as WEIGHTS array.
-         ! weights_(1)=wsc
+          weights_(1)=wsc
           weights_(2)=wscp
           weights_(3)=welec
           weights_(4)=wcorr
@@ -610,9 +610,9 @@
              .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0 &
              .or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then
 #endif
-!            write(iout,*),"just befor eelec call"
+!            print *,"just befor eelec call"
             call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
-!         write (iout,*) "ELEC calc"
+!            print *, "ELEC calc"
          else
             ees=0.0d0
             evdw1=0.0d0
@@ -831,6 +831,7 @@
        eespp=0.0d0
       endif
 !      write(iout,*) ecorr_nucl,"ecorr_nucl",nres_molec(2)
+!      print *,"before ecatcat"
       if (nfgtasks.gt.1) then
       if (fg_rank.eq.0) then
       call ecatcat(ecationcation)
@@ -2727,13 +2728,18 @@
 ! to calculate the el-loc multibody terms of various order.
 !
 !AL el      mu=0.0d0
+   
 #ifdef PARMAT
       do i=ivec_start+2,ivec_end+2
 #else
       do i=3,nres+1
 #endif
         if (i.gt. nnt+2 .and. i.lt.nct+2) then
+          if (itype(i-2,1).eq.0) then 
+          iti = nloctyp
+          else
           iti = itype2loc(itype(i-2,1))
+          endif
         else
           iti=nloctyp
         endif
@@ -3025,17 +3031,17 @@
 !d        write (iout,*) 'mu2',mu2(:,i-2)
         if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) &
         then  
-        call matmat2(CC(1,1,i-2),Ugder(1,1,i-2),CUgder(1,1,i-2))
+        call matmat2(CC(1,1,i-1),Ugder(1,1,i-2),CUgder(1,1,i-2))
         call matmat2(DD(1,1,i-2),Ugder(1,1,i-2),DUgder(1,1,i-2))
         call matmat2(Dtilde(1,1,i-2),Ug2der(1,1,i-2),DtUg2der(1,1,i-2))
-        call matvec2(Ctilde(1,1,i-2),obrot_der(1,i-2),Ctobrder(1,i-2))
+        call matvec2(Ctilde(1,1,i-1),obrot_der(1,i-2),Ctobrder(1,i-2))
         call matvec2(Dtilde(1,1,i-2),obrot2_der(1,i-2),Dtobr2der(1,i-2))
 ! Vectors and matrices dependent on a single virtual-bond dihedral.
         call matvec2(DD(1,1,i-2),b1tilde(1,iti1),auxvec(1))
         call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2)) 
         call matvec2(Ug2der(1,1,i-2),auxvec(1),Ug2Db1tder(1,i-2)) 
-        call matvec2(CC(1,1,i-2),Ub2(1,i-2),CUgb2(1,i-2))
-        call matvec2(CC(1,1,i-2),Ub2der(1,i-2),CUgb2der(1,i-2))
+        call matvec2(CC(1,1,i-1),Ub2(1,i-2),CUgb2(1,i-2))
+        call matvec2(CC(1,1,i-1),Ub2der(1,i-2),CUgb2der(1,i-2))
         call matmat2(EUg(1,1,i-2),CC(1,1,iti1),EUgC(1,1,i-2))
         call matmat2(EUgder(1,1,i-2),CC(1,1,iti1),EUgCder(1,1,i-2))
         call matmat2(EUg(1,1,i-2),DD(1,1,iti1),EUgD(1,1,i-2))
@@ -4338,7 +4344,9 @@
           +a23*gmuij1(2)&
           +a32*gmuij1(3)&
           +a33*gmuij1(4))&
-         *fac_shield(i)*fac_shield(j)
+         *fac_shield(i)*fac_shield(j)&
+                    *sss_ele_cut
+
 !c         write(iout,*) "derivative over thatai"
 !c         write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
 !c     &   a33*gmuij1(4) 
@@ -4354,7 +4362,9 @@
           +a33*gmuij2(4)
          gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+&
            geel_loc_ij*wel_loc&
-         *fac_shield(i)*fac_shield(j)
+         *fac_shield(i)*fac_shield(j)&
+                    *sss_ele_cut
+
 
 !c  Derivative over j residue
          geel_loc_ji=a22*gmuji1(1)&
@@ -4367,7 +4377,9 @@
 
         gloc(nphi+j,icg)=gloc(nphi+j,icg)+&
            geel_loc_ji*wel_loc&
-         *fac_shield(i)*fac_shield(j)
+         *fac_shield(i)*fac_shield(j)&
+                    *sss_ele_cut
+
 
          geel_loc_ji=&
           +a22*gmuji2(1)&
@@ -4379,7 +4391,8 @@
 !c     &   a33*gmuji2(4)
          gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+&
            geel_loc_ji*wel_loc&
-         *fac_shield(i)*fac_shield(j)
+         *fac_shield(i)*fac_shield(j)&
+                    *sss_ele_cut
 #endif
 
 !          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
@@ -11194,6 +11207,7 @@
 #ifdef TIMING
       time01=MPI_Wtime()
 #endif
+!#define DEBUG
 #ifdef DEBUG
       write (iout,*) "sum_gradient gvdwc, gvdwx"
       do i=1,nres
@@ -11227,18 +11241,18 @@
 !        write (iout,'(i5,3f10.5,2x,f10.5)') 
 !     &  i,(gcorr4_turn(j,i),j=1,3),gel_loc_turn4(i)
 !      enddo
-      write (iout,*) "gvdwc gvdwc_scp gvdwc_scpp"
-      do i=1,nres
-        write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
-         i,(gvdwc(j,i),j=1,3),(gvdwc_scp(j,i),j=1,3),&
-         (gvdwc_scpp(j,i),j=1,3)
-      enddo
-      write (iout,*) "gelc_long gvdwpp gel_loc_long"
-      do i=1,nres
-        write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
-         i,(gelc_long(j,i),j=1,3),(gvdwpp(j,i),j=1,3),&
-         (gelc_loc_long(j,i),j=1,3)
-      enddo
+!      write (iout,*) "gvdwc gvdwc_scp gvdwc_scpp"
+!      do i=1,nres
+!        write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
+!         i,(gvdwc(j,i),j=1,3),(gvdwc_scp(j,i),j=1,3),&
+!         (gvdwc_scpp(j,i),j=1,3)
+!      enddo
+!      write (iout,*) "gelc_long gvdwpp gel_loc_long"
+!      do i=1,nres
+!        write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') &
+!         i,(gelc_long(j,i),j=1,3),(gvdwpp(j,i),j=1,3),&
+!         (gelc_loc_long(j,i),j=1,3)
+!      enddo
       call flush(iout)
 #endif
 #ifdef SPLITELE
@@ -12749,7 +12763,7 @@
 !      call intcartderiv
 !      call checkintcartgrad
       call zerograd
-      aincr=2.0D-5
+      aincr=1.0D-7
       write(iout,*) 'Calling CHECK_ECARTINT.',aincr
       nf=0
       icall=0
@@ -22654,7 +22668,8 @@
         rs,rthrp,rfourp,rsixp,reight,Irsixp,Ireight,Irtw,Irfourt,      &
         opt1,opt2,opt3,opt4,opt5,opt6,opt7,opt8,opt9,opt10,opt11,opt12,&
         opt13,opt14,opt15,opt16,opt17,opt18,opt19, &
-        Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip
+        Equad1,Equad2,dscmag,v1dpv2,dscmag3,constA,constB,Edip,&
+        ndiv,ndivi
         real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,&
         gg,r,EtotalCat,dEtotalCm,dEtotalCalp,dEvan1Cm,dEvan2Cm, &
         dEtotalpep,dEtotalcat_num,dEddci,dEtotalcm_num,dEtotalcalp_num, &
@@ -22668,17 +22683,6 @@
         ecation_prot=0.0d0
 ! first lets calculate interaction with peptide groups
         if (nres_molec(5).eq.0) return
-         wconst=78
-        wdip =1.092777950857032D2
-        wdip=wdip/wconst
-        wmodquad=-2.174122713004870D4
-        wmodquad=wmodquad/wconst
-        wquad1 = 3.901232068562804D1
-        wquad1=wquad1/wconst
-        wquad2 = 3
-        wquad2=wquad2/wconst
-        wvan1 = 0.1
-        wvan2 = 6
         itmp=0
         do i=1,4
         itmp=itmp+nres_molec(i)
@@ -22698,6 +22702,23 @@
           if (zi.lt.0) zi=zi+boxzsize
 
          do j=itmp+1,itmp+nres_molec(5)
+!           print *,"WTF",itmp,j,i
+! all parameters were for Ca2+ to approximate single charge divide by two
+         ndiv=1.0
+         if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+         wconst=78*ndiv
+        wdip =1.092777950857032D2
+        wdip=wdip/wconst
+        wmodquad=-2.174122713004870D4
+        wmodquad=wmodquad/wconst
+        wquad1 = 3.901232068562804D1
+        wquad1=wquad1/wconst
+        wquad2 = 3
+        wquad2=wquad2/wconst
+        wvan1 = 0.1
+        wvan2 = 6
+!        itmp=0
+
            xj=c(1,j)
            yj=c(2,j)
            zj=c(3,j)
@@ -22771,6 +22792,7 @@
         E2 = -wquad1*Irthrp*wquad2+wvan1*(wvan2**12*Irtwelv-&
              2*wvan2**6*Irsistp)
         ecation_prot = ecation_prot+E1+E2
+!        print *,"ecatprot",i,j,ecation_prot,rcpm
         dE1dr = -2*costhet*wdip*Irthrp-& 
          (4*wmodquad*Irfiftp+3*wquad1*Irfourp)*sin2thet
         dE2dr = 3*wquad1*wquad2*Irfourp-     &
@@ -22810,6 +22832,9 @@
           enddo
            cm1mag=sqrt(cm1(1)**2+cm1(2)**2+cm1(3)**2)
          do j=itmp+1,itmp+nres_molec(5)
+         ndiv=1.0
+         if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
            xj=c(1,j)
            yj=c(2,j)
            zj=c(3,j)
@@ -22852,7 +22877,10 @@
        endif
 !       enddo
 !       enddo
-         if(itype(i,1).eq.15.or.itype(i,1).eq.16) then
+! 15- Glu 16-Asp
+         if((itype(i,1).eq.15.or.itype(i,1).eq.16).or.&
+         ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.&
+         (itype(i,1).eq.25))) then
             if(itype(i,1).eq.16) then
             inum=1
             else
@@ -22897,7 +22925,15 @@
 
 !  The weights of the energy function calculated from
 !The quantum mechanical GAMESS simulations of calcium with ASP/GLU
-        wh2o=78
+          if ((itype(i,1).eq.27).or.(itype(i,1).eq.26).or.(itype(i,1).eq.25)) then
+            ndivi=0.5
+          else
+            ndivi=1.0
+          endif
+         ndiv=1.0
+         if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
+        wh2o=78*ndivi*ndiv
         wc = vcatprm(1)
         wc=wc/wh2o
         wdip =vcatprm(2)
@@ -23066,7 +23102,10 @@
         v1dpv2 = v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)
 !  The weights of the energy function calculated from
 !The quantum mechanical GAMESS simulations of ASN/GLN with calcium
-        wh2o=78
+         ndiv=1.0
+         if ((itype(j,5).eq.1).or.(itype(j,5).eq.3)) ndiv=2.0
+
+        wh2o=78*ndiv
         wdip =vcatprm(2)
         wdip=wdip/wh2o
         wquad1 =vcatprm(3)

diff --git a/source/unres/io.F90 b/source/unres/io.F90
index abce524..4f526c4 100644 (file)
--- a/source/unres/io.F90
+++ b/source/unres/io.F90
@@ -955,6 +955,8 @@
   34    continue
 !        print *,'Begin reading pdb data'
         call readpdb
+!        call int_from_cart1(.true.)
+
 !        print *,'Finished reading pdb data'
         if(me.eq.king.or..not.out1file) &
          write (iout,'(a,i3,a,i3)')'nsup=',nsup,&
@@ -1054,7 +1056,7 @@
        nres=0
        do i=1,5
         nres=nres+nres_molec(i)
-        print *,nres_molec(i)
+        print *,"nres_molec",nres,nres_molec(i)
        enddo
        
 ! Assign initial virtual bond lengths
@@ -1069,14 +1071,18 @@
          enddo
 !        print *,nres_molec(i)
         endif
-        if(.not.allocated(vbld)) allocate(vbld(2*nres))
+        print *,nres,"nres"
+        if(.not.allocated(vbld)) then
+           print *, "I DO ENTER" 
+           allocate(vbld(2*nres))
+        endif
         if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres))
         do i=2,nres
           vbld(i)=vbl
           vbld_inv(i)=vblinv
         enddo
         do i=2,nres-1
-          print *, "molnum",molnum(i),itype(i,molnum(i)) 
+          print *, "molnum",molnum(i),itype(i,molnum(i)),nres,i 
           vbld(i+nres)=dsc(iabs(itype(i,molnum(i))))
           vbld_inv(i+nres)=dsc_inv(iabs(itype(i,molnum(i))))
 !          write (iout,*) "i",i," itype",itype(i,1),
@@ -1553,6 +1559,10 @@
           enddo
       endif
       if (i2ndstr.gt.0) call secstrp2dihc
+      if (indpdb.gt.0) then 
+          write(iout,*) "WCHODZE TU!!"
+          call int_from_cart1(.true.)
+      endif
 !      call geom_to_var(nvar,x)
 !      call etotal(energia(0))
 !      call enerprint(energia(0))

diff --git a/source/unres/io_config.F90 b/source/unres/io_config.F90
index 5907752..b9a9844 100644 (file)
--- a/source/unres/io_config.F90
+++ b/source/unres/io_config.F90
@@ -1573,6 +1573,7 @@
         itype2loc(-i)=-itype2loc(i)
       enddo
 #else
+      allocate(iloctyp(-nloctyp:nloctyp))
       iloctyp(0)=10
       iloctyp(1)=9
       iloctyp(2)=20
@@ -1825,7 +1826,7 @@
       allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
       allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
       allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
-
+      allocate(b(13,-nloctyp-1:nloctyp+1))
       if (lprint) &
        write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
       do i=0,nloctyp-1
@@ -1841,35 +1842,35 @@
         b(4,-i)=-b(4,i)
         b(5,-i)=-b(5,i)
         endif
-c        B1(1,i)  = b(3)
-c        B1(2,i)  = b(5)
-c        B1(1,-i) = b(3)
-c        B1(2,-i) = -b(5)
-c        b1(1,i)=0.0d0
-c        b1(2,i)=0.0d0
-c        B1tilde(1,i) = b(3)
-c        B1tilde(2,i) =-b(5)
-c        B1tilde(1,-i) =-b(3)
-c        B1tilde(2,-i) =b(5)
-c        b1tilde(1,i)=0.0d0
-c        b1tilde(2,i)=0.0d0
-c        B2(1,i)  = b(2)
-c        B2(2,i)  = b(4)
-c        B2(1,-i)  =b(2)
-c        B2(2,-i)  =-b(4)
-cc        B1tilde(1,i) = b(3,i)
-cc        B1tilde(2,i) =-b(5,i)
-C        B1tilde(1,-i) =-b(3,i)
-C        B1tilde(2,-i) =b(5,i)
-cc        b1tilde(1,i)=0.0d0
-cc        b1tilde(2,i)=0.0d0
-cc        B2(1,i)  = b(2,i)
-cc        B2(2,i)  = b(4,i)
-C        B2(1,-i)  =b(2,i)
-C        B2(2,-i)  =-b(4,i)
-
-c        b2(1,i)=0.0d0
-c        b2(2,i)=0.0d0
+!c        B1(1,i)  = b(3)
+!c        B1(2,i)  = b(5)
+!c        B1(1,-i) = b(3)
+!c        B1(2,-i) = -b(5)
+!c        b1(1,i)=0.0d0
+!c        b1(2,i)=0.0d0
+!c        B1tilde(1,i) = b(3)
+!c!        B1tilde(2,i) =-b(5)
+!c!        B1tilde(1,-i) =-b(3)
+!c!        B1tilde(2,-i) =b(5)
+!c!        b1tilde(1,i)=0.0d0
+!c        b1tilde(2,i)=0.0d0
+!c        B2(1,i)  = b(2)
+!c        B2(2,i)  = b(4)
+!c        B2(1,-i)  =b(2)
+!c        B2(2,-i)  =-b(4)
+!cc        B1tilde(1,i) = b(3,i)
+!cc        B1tilde(2,i) =-b(5,i)
+!c        B1tilde(1,-i) =-b(3,i)
+!c        B1tilde(2,-i) =b(5,i)
+!cc        b1tilde(1,i)=0.0d0
+!cc        b1tilde(2,i)=0.0d0
+!cc        B2(1,i)  = b(2,i)
+!cc        B2(2,i)  = b(4,i)
+!c        B2(1,-i)  =b(2,i)
+!c        B2(2,-i)  =-b(4,i)
+
+!c        b2(1,i)=0.0d0
+!c        b2(2,i)=0.0d0
         CCold(1,1,i)= b(7,i)
         CCold(2,2,i)=-b(7,i)
         CCold(2,1,i)= b(9,i)
@@ -1878,27 +1879,27 @@ c        b2(2,i)=0.0d0
         CCold(2,2,-i)=-b(7,i)
         CCold(2,1,-i)=-b(9,i)
         CCold(1,2,-i)=-b(9,i)
-c        CC(1,1,i)=0.0d0
-c        CC(2,2,i)=0.0d0
-c        CC(2,1,i)=0.0d0
-c        CC(1,2,i)=0.0d0
-c        Ctilde(1,1,i)= CCold(1,1,i)
-c        Ctilde(1,2,i)= CCold(1,2,i)
-c        Ctilde(2,1,i)=-CCold(2,1,i)
-c        Ctilde(2,2,i)=-CCold(2,2,i)
-c        CC(1,1,i)=0.0d0
-c        CC(2,2,i)=0.0d0
-c        CC(2,1,i)=0.0d0
-c        CC(1,2,i)=0.0d0
-c        Ctilde(1,1,i)= CCold(1,1,i)
-c        Ctilde(1,2,i)= CCold(1,2,i)
-c        Ctilde(2,1,i)=-CCold(2,1,i)
-c        Ctilde(2,2,i)=-CCold(2,2,i)
-
-c        Ctilde(1,1,i)=0.0d0
-c        Ctilde(1,2,i)=0.0d0
-c        Ctilde(2,1,i)=0.0d0
-c        Ctilde(2,2,i)=0.0d0
+!c        CC(1,1,i)=0.0d0
+!c        CC(2,2,i)=0.0d0
+!c        CC(2,1,i)=0.0d0
+!c        CC(1,2,i)=0.0d0
+!c        Ctilde(1,1,i)= CCold(1,1,i)
+!c        Ctilde(1,2,i)= CCold(1,2,i)
+!c        Ctilde(2,1,i)=-CCold(2,1,i)
+!c        Ctilde(2,2,i)=-CCold(2,2,i)
+!c        CC(1,1,i)=0.0d0
+!c        CC(2,2,i)=0.0d0
+!c        CC(2,1,i)=0.0d0
+!c        CC(1,2,i)=0.0d0
+!c        Ctilde(1,1,i)= CCold(1,1,i)
+!c        Ctilde(1,2,i)= CCold(1,2,i)
+!c        Ctilde(2,1,i)=-CCold(2,1,i)
+!c        Ctilde(2,2,i)=-CCold(2,2,i)
+
+!c        Ctilde(1,1,i)=0.0d0
+!c        Ctilde(1,2,i)=0.0d0
+!c        Ctilde(2,1,i)=0.0d0
+!c        Ctilde(2,2,i)=0.0d0
         DDold(1,1,i)= b(6,i)
         DDold(2,2,i)=-b(6,i)
         DDold(2,1,i)= b(8,i)
@@ -1907,19 +1908,19 @@ c        Ctilde(2,2,i)=0.0d0
         DDold(2,2,-i)=-b(6,i)
         DDold(2,1,-i)=-b(8,i)
         DDold(1,2,-i)=-b(8,i)
-c        DD(1,1,i)=0.0d0
-c        DD(2,2,i)=0.0d0
-c        DD(2,1,i)=0.0d0
-c        DD(1,2,i)=0.0d0
-c        Dtilde(1,1,i)= DD(1,1,i)
-c        Dtilde(1,2,i)= DD(1,2,i)
-c        Dtilde(2,1,i)=-DD(2,1,i)
-c        Dtilde(2,2,i)=-DD(2,2,i)
-
-c        Dtilde(1,1,i)=0.0d0
-c        Dtilde(1,2,i)=0.0d0
-c        Dtilde(2,1,i)=0.0d0
-c        Dtilde(2,2,i)=0.0d0
+!c        DD(1,1,i)=0.0d0
+!c        DD(2,2,i)=0.0d0
+!c        DD(2,1,i)=0.0d0
+!c        DD(1,2,i)=0.0d0
+!c        Dtilde(1,1,i)= DD(1,1,i)
+!c        Dtilde(1,2,i)= DD(1,2,i)
+!c        Dtilde(2,1,i)=-DD(2,1,i)
+!c        Dtilde(2,2,i)=-DD(2,2,i)
+
+!c        Dtilde(1,1,i)=0.0d0
+!c        Dtilde(1,2,i)=0.0d0
+!c        Dtilde(2,1,i)=0.0d0
+!c        Dtilde(2,2,i)=0.0d0
         EEold(1,1,i)= b(10,i)+b(11,i)
         EEold(2,2,i)=-b(10,i)+b(11,i)
         EEold(2,1,i)= b(12,i)-b(13,i)
@@ -1930,11 +1931,11 @@ c        Dtilde(2,2,i)=0.0d0
         EEold(1,2,-i)=-b(12,i)-b(13,i)
         write(iout,*) "TU DOCHODZE"
         print *,"JESTEM"
-c        ee(1,1,i)=1.0d0
-c        ee(2,2,i)=1.0d0
-c        ee(2,1,i)=0.0d0
-c        ee(1,2,i)=0.0d0
-c        ee(2,1,i)=ee(1,2,i)
+!c        ee(1,1,i)=1.0d0
+!c        ee(2,2,i)=1.0d0
+!c        ee(2,1,i)=0.0d0
+!c        ee(1,2,i)=0.0d0
+!c        ee(2,1,i)=ee(1,2,i)
       enddo
       if (lprint) then
       write (iout,*)
@@ -3407,6 +3408,8 @@ c        ee(2,1,i)=ee(1,2,i)
         if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
 ! Fish out the ATOM cards.
 !        write(iout,*) 'card',card(1:20)
+!        print *,"ATU ",card(1:6), CARD(3:6)
+!        print *,card
         if (index(card(1:4),'ATOM').gt.0) then  
           read (card(12:16),*) atom
 !          write (iout,*) "! ",atom," !",ires
@@ -3428,9 +3431,9 @@ c        ee(2,1,i)=ee(1,2,i)
                 enddo
               else
                 call sccenter(ires_old,iii,sccor)
-              endif
+              endif !unres_pdb
               iii=0
-            endif
+            endif !ind_pdb
 ! Start new residue.
             if (res.eq.'Cl-' .or. res.eq.'Na+') then
               ires=ires_old
@@ -3442,7 +3445,7 @@ c        ee(2,1,i)=ee(1,2,i)
                 ishift=ishift-1
                 itype(1,1)=ntyp1
                 nres_molec(molecule)=nres_molec(molecule)+1
-              endif
+              endif ! Gly
               ires=ires-ishift+ishift1
               ires_old=ires
 !              write (iout,*) "ishift",ishift," ires",ires,&
@@ -3461,7 +3464,7 @@ c        ee(2,1,i)=ee(1,2,i)
               ishift=ishift-(ires-ishift+ishift1-ires_old-1)
               ires=ires-ishift+ishift1
               ires_old=ires
-            endif 
+            endif ! Na Cl
 !            print *,'atom',ires,atom
             if (res.eq.'ACE' .or. res.eq.'NHE') then
               itype(ires,1)=10
@@ -3485,11 +3488,11 @@ c        ee(2,1,i)=ee(1,2,i)
 !              print *,"ires=",ires,istype(ires)
 !            endif
 
-            endif
-            endif
+            endif ! atom.eq.CA
+            endif !ACE
           else
             ires=ires-ishift+ishift1
-          endif
+          endif !ires_old
 !          write (iout,*) "ires_old",ires_old," ires",ires
           if (card(27:27).eq."A" .or. card(27:27).eq."B") then
 !            ishift1=ishift1+1
@@ -3574,6 +3577,7 @@ c        ee(2,1,i)=ee(1,2,i)
             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
               endif
           endif
+!         print *,"IONS",ions,card(1:6)
         else if ((ions).and.(card(1:6).eq.'HETATM')) then
        if (firstion.eq.0) then 
        firstion=1
@@ -3583,10 +3587,10 @@ c        ee(2,1,i)=ee(1,2,i)
          enddo
        else
           call sccenter(ires,iii,sccor)
-       endif
-       endif
+       endif ! unres_pdb
+       endif !firstion
           read (card(12:16),*) atom
-          print *,"HETATOM", atom
+!          print *,"HETATOM", atom
           read (card(18:20),'(a3)') res
           if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
           (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG')           &
@@ -3600,7 +3604,7 @@ c        ee(2,1,i)=ee(1,2,i)
            res=res(2:3)//' '
            itype(ires,molecule)=rescode(ires,res,0,molecule)
            read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
-          endif
+          endif! NA
         endif !atom
       enddo
    10 write (iout,'(a,i5)') ' Number of residues found: ',ires
@@ -3674,7 +3678,8 @@ c        ee(2,1,i)=ee(1,2,i)
 !              e2(3)=0.0d0
 !            endif
             do j=1,3
-              c(j,i)=c(j,i+1)-1.9d0*e2(j)
+!              c(j,i)=c(j,i+1)-1.9d0*e2(j)
+             c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
             enddo
 !          else !unres_pdb
            do j=1,3
@@ -3789,7 +3794,8 @@ c        ee(2,1,i)=ee(1,2,i)
             e2(3)=0.0d0
           endif
           do j=1,3
-            c(j,1)=c(j,2)-1.9d0*e2(j)
+!            c(j,1)=c(j,2)-1.9d0*e2(j)
+             c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
           enddo
         else
         do j=1,3

diff --git a/source/unres/unres.F90 b/source/unres/unres.F90
index bc9da96..4a8e7a6 100644 (file)
--- a/source/unres/unres.F90
+++ b/source/unres/unres.F90
@@ -818,6 +818,7 @@
       call cartprint
       call intout
       icall=1
+      write (iout,*) "before etotal"
       call etotal(energy_(0))
       etot = energy_(0)
       call enerprint(energy_(0))