dzialajacy gradient

Conflicts:

	source/unres/src_MD-M/energy_p_new_barrier.F

diff --git a/source/unres/src_MD-M/CMakeLists.txt b/source/unres/src_MD-M/CMakeLists.txt
index cfdc654..3884bc1 100644 (file)
--- a/source/unres/src_MD-M/CMakeLists.txt
+++ b/source/unres/src_MD-M/CMakeLists.txt
@@ -203,7 +203,7 @@ set_property(SOURCE ${UNRES_MDM_SRC3} PROPERTY COMPILE_FLAGS ${FFLAGS3} )
 #=========================================
 if(UNRES_MD_FF STREQUAL "GAB" )
   # set preprocesor flags   
-  set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC  -DSCCORPDB" )
+  set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC  -DSCCORPDB -DNEWCORR" )
 
 #=========================================
 #  Settings for E0LL2Y force field

diff --git a/source/unres/src_MD-M/COMMON.CONTACTS b/source/unres/src_MD-M/COMMON.CONTACTS
index 4b81714..45c578b 100644 (file)
--- a/source/unres/src_MD-M/COMMON.CONTACTS
+++ b/source/unres/src_MD-M/COMMON.CONTACTS
@@ -28,7 +28,8 @@ C 10/30/99 Added other pre-computed vectors and matrices needed
 C          to calculate three - six-order el-loc correlation terms
       double precision Ug,Ugder,Ug2,Ug2der,obrot,obrot2,obrot_der,
      &  obrot2_der,Ub2,Ub2der,mu,muder,EUg,EUgder,CUg,CUgder,gmu,gUb2
-     &  DUg,DUgder,DtUg2,DtUg2der,Ctobr,Ctobrder,Dtobr2,Dtobr2der
+     &  DUg,DUgder,DtUg2,DtUg2der,Ctobr,Ctobrder,Dtobr2,Dtobr2der,
+     &  gtEug
       common /rotat/ Ug(2,2,maxres),Ugder(2,2,maxres),Ug2(2,2,maxres),
      &  Ug2der(2,2,maxres),obrot(2,maxres),obrot2(2,maxres),
      &  obrot_der(2,maxres),obrot2_der(2,maxres)
@@ -40,7 +41,7 @@ C amino-acid residue.
      &  Dtobr2(2,maxres),Dtobr2der(2,maxres),
      &  EUg(2,2,maxres),EUgder(2,2,maxres),CUg(2,2,maxres),
      &  CUgder(2,2,maxres),DUg(2,2,maxres),Dugder(2,2,maxres),
-     &  DtUg2(2,2,maxres),DtUg2der(2,2,maxres)
+     &  DtUg2(2,2,maxres),DtUg2der(2,2,maxres),gtEUg(2,2,maxres)
 C This common block contains vectors and matrices dependent on two
 C consecutive amino-acid residues.
       double precision Ug2Db1t,Ug2Db1tder,CUgb2,CUgb2der,EUgC,

diff --git a/source/unres/src_MD-M/COMMON.TORSION b/source/unres/src_MD-M/COMMON.TORSION
index 95472be..3f1b981 100644 (file)
--- a/source/unres/src_MD-M/COMMON.TORSION
+++ b/source/unres/src_MD-M/COMMON.TORSION
@@ -25,16 +25,17 @@ C 6/23/01 - constants for double torsionals
 C 9/18/99 - added Fourier coeffficients of the expansion of local energy 
 C           surface
       double precision b1,b2,cc,dd,ee,ctilde,dtilde,b2tilde,b1tilde,
-     &bnew1,bnew2,eenew,ggb1,ggb2
+     &bnew1,bnew2,eenew,gtb1,gtb2,eeold,gtee
       integer nloctyp
       common/fourier/ b1(2,maxres),b2(2,maxres),
      &    bnew1(3,2,-maxtor:maxtor),bnew2(3,2,-maxtor:maxtor),
      &    cc(2,2,-maxtor:maxtor),
-     &    dd(2,2,-maxtor:maxtor),ee(2,2,-maxtor:maxtor),
+     &    dd(2,2,-maxtor:maxtor),eeold(2,2,-maxtor:maxtor),
      &    eenew(2,-maxtor:maxtor),
+     &    ee(2,2,maxres),
      &    ctilde(2,2,-maxtor:maxtor),
      &    dtilde(2,2,-maxtor:maxtor),b1tilde(2,maxres),
      &    b2tilde(2,maxres),
-     &    gtb1(2,maxres),gtb2(2,maxres),
+     &    gtb1(2,maxres),gtb2(2,maxres),gtEE(2,2,maxres),
      &    nloctyp
      

diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index c13d341..ab7d468 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -2297,7 +2297,14 @@ c        endif
         gtb1(2,i-2)=0.0
         b2(2,i-2)=bnew2(1,2,iti)
         gtb2(2,i-2)=0.0
-c        EE(1,1,iti)=0.0d0
+        EE(1,1,i-2)=eenew(1,iti)*dcos(theta(i-1))
+        EE(1,2,i-2)=eeold(1,2,iti)
+        EE(2,1,i-2)=eeold(2,1,iti)
+        EE(2,2,i-2)=eeold(2,2,iti)
+        gtEE(1,1,i-2)=-eenew(1,iti)*dsin(theta(i-1))
+        gtEE(1,2,i-2)=0.0d0
+        gtEE(2,2,i-2)=0.0d0
+        gtEE(2,1,i-2)=0.0d0
 c        EE(2,2,iti)=0.0d0
 c        EE(1,2,iti)=0.5d0*eenew(1,iti)
 c        EE(2,1,iti)=0.5d0*eenew(1,iti)
@@ -2383,7 +2390,6 @@ c        b2(2,iti)=bnew2(1,2,iti)*sin(alpha(i))*sin(beta(i))
           Ug2der(2,2,i-2)=0.0d0
         endif
 c        if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then
-#ifndef NEWCORR
         if (i.gt. nnt+2 .and. i.lt.nct+2) then
           iti = itortyp(itype(i-2))
         else
@@ -2395,7 +2401,6 @@ c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
         else
           iti1=ntortyp+1
         endif
-#endif
 cd        write (iout,*) '*******i',i,' iti1',iti
 cd        write (iout,*) 'b1',b1(:,iti)
 cd        write (iout,*) 'b2',b2(:,iti)
@@ -2407,7 +2412,13 @@ c        if (i .gt. iatel_s+2) then
           call matvec2(Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2))
 c          write (iout,*) Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2),"chuj"
 #endif
-          call matmat2(EE(1,1,iti),Ug(1,1,i-2),EUg(1,1,i-2))
+c          write(iout,*) "co jest kurwa", iti, EE(1,1,iti),EE(2,1,iti),
+c     &    EE(1,2,iti),EE(2,2,iti)
+          call matmat2(EE(1,1,i-2),Ug(1,1,i-2),EUg(1,1,i-2))
+          call matmat2(gtEE(1,1,i-2),Ug(1,1,i-2),gtEUg(1,1,i-2))
+c          write(iout,*) "Macierz EUG",
+c     &    eug(1,1,i-2),eug(1,2,i-2),eug(2,1,i-2),
+c     &    eug(2,2,i-2)
           if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) 
      &    then
           call matmat2(CC(1,1,iti),Ug(1,1,i-2),CUg(1,1,i-2))
@@ -2430,7 +2441,7 @@ c          write (iout,*) Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2),"chuj"
           enddo
         endif
         call matvec2(Ugder(1,1,i-2),b2(1,i-2),Ub2der(1,i-2))
-        call matmat2(EE(1,1,iti),Ugder(1,1,i-2),EUgder(1,1,i-2))
+        call matmat2(EE(1,1,i-2),Ugder(1,1,i-2),EUgder(1,1,i-2))
         do k=1,2
           muder(k,i-2)=Ub2der(k,i-2)
         enddo
@@ -2866,8 +2877,7 @@ C
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
         num_conti=0
-c TU ZLE
-c        call eelecij(i,i+2,ees,evdw1,eel_loc)
+        call eelecij(i,i+2,ees,evdw1,eel_loc)
         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
         num_cont_hb(i)=num_conti
       enddo
@@ -2885,8 +2895,8 @@ c        call eelecij(i,i+2,ees,evdw1,eel_loc)
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
         num_conti=num_cont_hb(i)
-c TU ZLE
-c        call eelecij(i,i+3,ees,evdw1,eel_loc)
+c        write(iout,*) "JESTEM W PETLI"
+        call eelecij(i,i+3,ees,evdw1,eel_loc)
         if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) 
      &   call eturn4(i,eello_turn4)
         num_cont_hb(i)=num_conti
@@ -2894,8 +2904,8 @@ c        call eelecij(i,i+3,ees,evdw1,eel_loc)
 c
 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 c
-c      do i=iatel_s,iatel_e
-       do i=4,5
+      do i=iatel_s,iatel_e
+c       do i=7,7
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
@@ -2908,9 +2918,9 @@ c      do i=iatel_s,iatel_e
         zmedi=c(3,i)+0.5d0*dzi
 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         num_conti=num_cont_hb(i)
-c        do j=ielstart(i),ielend(i)
-         do j=8,9
-          write (iout,*) 'tu wchodze',i,j,itype(i),itype(j)
+        do j=ielstart(i),ielend(i)
+c         do j=13,13
+c          write (iout,*) 'tu wchodze',i,j,itype(i),itype(j)
           if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle
           call eelecij(i,j,ees,evdw1,eel_loc)
         enddo ! j
@@ -3640,7 +3650,9 @@ C Third- and fourth-order contributions from turns
       dimension ggg(3)
       double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2),
      &  e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2),
-     &  e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2)
+     &  e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),gpizda1(2,2),
+     &  gpizda2(2,2),auxgmat1(2,2),auxgmatt1(2,2),
+     &  auxgmat2(2,2),auxgmatt2(2,2)
       double precision agg(3,4),aggi(3,4),aggi1(3,4),
      &    aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2)
       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
@@ -3664,9 +3676,24 @@ C
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC   
 cd        call checkint_turn3(i,a_temp,eello_turn3_num)
         call matmat2(EUg(1,1,i+1),EUg(1,1,i+2),auxmat(1,1))
+c auxalary matices for theta gradient
+c auxalary matrix for i+1 and constant i+2
+        call matmat2(gtEUg(1,1,i+1),EUg(1,1,i+2),auxgmat1(1,1))
+c auxalary matrix for i+2 and constant i+1
+        call matmat2(EUg(1,1,i+1),gtEUg(1,1,i+2),auxgmat2(1,1))
         call transpose2(auxmat(1,1),auxmat1(1,1))
+        call transpose2(auxgmat1(1,1),auxgmatt1(1,1))
+        call transpose2(auxgmat2(1,1),auxgmatt2(1,1))
         call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1))
+        call matmat2(a_temp(1,1),auxgmatt1(1,1),gpizda1(1,1))
+        call matmat2(a_temp(1,1),auxgmatt2(1,1),gpizda2(1,1))
         eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2))
+C Derivatives in theta
+        gloc(nphi+i,icg)=gloc(nphi+i,icg)
+     &  +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3
+        gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)
+     &  +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3
+
         if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
      &          'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2))
 cd        write (2,*) 'i,',i,' j',j,'eello_turn3',
@@ -3740,7 +3767,11 @@ C Third- and fourth-order contributions from turns
       dimension ggg(3)
       double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2),
      &  e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2),
-     &  e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),auxgvec(2)
+     &  e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),auxgvec(2),
+     &  auxgEvec1(2),auxgEvec2(2),auxgEvec3(2),
+     &  gte1t(2,2),gte2t(2,2),gte3t(2,2),
+     &  gte1a(2,2),gtae3(2,2),gtae3e2(2,2), ae3gte2(2,2),
+     &  gtEpizda1(2,2),gtEpizda2(2,2),gtEpizda3(2,2)
       double precision agg(3,4),aggi(3,4),aggi1(3,4),
      &    aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2)
       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
@@ -3760,6 +3791,7 @@ C
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC   
 cd        call checkint_turn4(i,a_temp,eello_turn4_num)
 c        write (iout,*) "eturn4 i",i," j",j," j1",j1," j2",j2
+c        write(iout,*)"WCHODZE W PROGRAM"
         a_temp(1,1)=a22
         a_temp(1,2)=a23
         a_temp(2,1)=a32
@@ -3771,37 +3803,83 @@ c        write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3
         call transpose2(EUg(1,1,i+1),e1t(1,1))
         call transpose2(Eug(1,1,i+2),e2t(1,1))
         call transpose2(Eug(1,1,i+3),e3t(1,1))
+C Ematrix derivative in theta
+        call transpose2(gtEUg(1,1,i+1),gte1t(1,1))
+        call transpose2(gtEug(1,1,i+2),gte2t(1,1))
+        call transpose2(gtEug(1,1,i+3),gte3t(1,1))
         call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
+c       eta1 in derivative theta
+        call matmat2(gte1t(1,1),a_temp(1,1),gte1a(1,1))
         call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
+c       auxgvec is derivative of Ub2 so i+3 theta
         call matvec2(e1a(1,1),gUb2(1,i+3),auxgvec(1)) 
+c       auxalary matrix of E i+1
+        call matvec2(gte1a(1,1),Ub2(1,i+3),auxgEvec1(1))
 c        s1=0.0
 c        gs1=0.0    
         s1=scalar2(b1(1,i+2),auxvec(1))
-c        gs1=scalar2(gtb1(1,i+2),auxgvec(1))
+c derivative of theta i+2 with constant i+3
+        gs23=scalar2(gtb1(1,i+2),auxvec(1))
+c derivative of theta i+2 with constant i+2
+        gs32=scalar2(b1(1,i+2),auxgvec(1))
+c derivative of E matix in theta of i+1
+        gsE13=scalar2(b1(1,i+2),auxgEvec1(1))
+
         call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
+c       ea31 in derivative theta
+        call matmat2(a_temp(1,1),gte3t(1,1),gtae3(1,1))
         call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) 
+c auxilary matrix auxgvec of Ub2 with constant E matirx
         call matvec2(ae3(1,1),gUb2(1,i+2),auxgvec(1))
+c auxilary matrix auxgEvec1 of E matix with Ub2 constant
+        call matvec2(gtae3(1,1),Ub2(1,i+2),auxgEvec3(1))
+
 c        s2=0.0
 c        gs2=0.0
         s2=scalar2(b1(1,i+1),auxvec(1))
-c        gs2=scalar2(gtb1(1,i+1),auxgvec(1))
-c        write(iout,*) gs1,gs2,'i=',i,auxgvec(1),gUb2(1,i+2),ggb1(1,i+2),
-c     &  ggb1(1,i+1)
+c derivative of theta i+1 with constant i+3
+        gs13=scalar2(gtb1(1,i+1),auxvec(1))
+c derivative of theta i+2 with constant i+1
+        gs21=scalar2(b1(1,i+1),auxgvec(1))
+c derivative of theta i+3 with constant i+1
+        gsE31=scalar2(b1(1,i+1),auxgEvec3(1))
+c        write(iout,*) gs1,gs2,'i=',i,auxgvec(1),gUb2(1,i+2),gtb1(1,i+2),
+c     &  gtb1(1,i+1)
         call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
+c two derivatives over diffetent matrices
+c gtae3e2 is derivative over i+3
+        call matmat2(gtae3(1,1),e2t(1,1),gtae3e2(1,1))
+c ae3gte2 is derivative over i+2
+        call matmat2(ae3(1,1),gte2t(1,1),ae3gte2(1,1))
         call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
+c three possible derivative over theta E matices
+c i+1
+        call matmat2(ae3e2(1,1),gte1t(1,1),gtEpizda1(1,1))
+c i+2
+        call matmat2(ae3gte2(1,1),e1t(1,1),gtEpizda2(1,1))
+c i+3
+        call matmat2(gtae3e2(1,1),e1t(1,1),gtEpizda3(1,1))
         s3=0.5d0*(pizda(1,1)+pizda(2,2))
+
+        gsEE1=0.5d0*(gtEpizda1(1,1)+gtEpizda1(2,2))
+        gsEE2=0.5d0*(gtEpizda2(1,1)+gtEpizda2(2,2))
+        gsEE3=0.5d0*(gtEpizda3(1,1)+gtEpizda3(2,2))
+
         eello_turn4=eello_turn4-(s1+s2+s3)
 #ifdef NEWCORR
-c        geel_loc_ij=-(gs1+gs2)
-c         gloc(nphi+i,icg)=gloc(nphi+i,icg)-
-c     &   gs1
+        gloc(nphi+i,icg)=gloc(nphi+i,icg)
+     &                  -(gs13+gsE13+gsEE1)*wturn4
+        gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg)
+     &                    -(gs23+gs21+gsEE2)*wturn4
+        gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg)
+     &                    -(gs32+gsE31+gsEE3)*wturn4
 c         gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)-
 c     &   gs2
 #endif
         if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
      &      'eturn4',i,j,-(s1+s2+s3)
-cd        write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
-cd     &    ' eello_turn4_num',8*eello_turn4_num
+c        write (iout,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
+c     &    ' eello_turn4_num',8*eello_turn4_num
 C Derivatives in gamma(i)
         call transpose2(EUgder(1,1,i+1),e1tder(1,1))
         call matmat2(e1tder(1,1),a_temp(1,1),auxmat(1,1))
@@ -4544,7 +4622,7 @@ C Derivatives of the "mean" values in gamma1 and gamma2.
      &      'ebend',i,ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
-        gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)
+        gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg)
       enddo
 C Ufff.... We've done all this!!! 
       return
@@ -4850,7 +4928,7 @@ c        lprn1=.false.
         etheta=etheta+ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
-        gloc(nphi+i-2,icg)=wang*dethetai
+        gloc(nphi+i-2,icg)=wang*dethetai+gloc(nphi+i-2,icg)
       enddo
       return
       end

diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F
index 9db9b60..03efdb2 100644 (file)
--- a/source/unres/src_MD-M/parmread.F
+++ b/source/unres/src_MD-M/parmread.F
@@ -991,21 +991,22 @@ 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
-        EE(1,1,i)= b(10)+b(11)
-        EE(2,2,i)=-b(10)+b(11)
-        EE(2,1,i)= b(12)-b(13)
-        EE(1,2,i)= b(12)+b(13)
-        EE(1,1,-i)= b(10)+b(11)
-        EE(2,2,-i)=-b(10)+b(11)
-        EE(2,1,-i)=-b(12)+b(13)
-        EE(1,2,-i)=-b(12)-b(13)
-
+        EEold(1,1,i)= b(10)+b(11)
+        EEold(2,2,i)=-b(10)+b(11)
+        EEold(2,1,i)= b(12)-b(13)
+        EEold(1,2,i)= b(12)+b(13)
+        EEold(1,1,-i)= b(10)+b(11)
+        EEold(2,2,-i)=-b(10)+b(11)
+        EEold(2,1,-i)=-b(12)+b(13)
+        EEold(1,2,-i)=-b(12)-b(13)
+        write(iout,*) "TU DOCHODZE"
 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
+c      lprint=.true.
       if (lprint) then
       do i=1,nloctyp
         write (iout,*) 'Type',i
@@ -1023,10 +1024,11 @@ c        ee(2,1,i)=ee(1,2,i)
         enddo
         write(iout,*) 'EE'
         do j=1,2
-          write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
+          write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
         enddo
       enddo
       endif
+c      lprint=.false.
 
 C 
 C Read electrostatic-interaction parameters