working microcanonical for CA2+ protein

diff --git a/PARAM/cation_param.parm b/PARAM/cation_param.parm
new file mode 100644 (file)
index 0000000..c6bf926
--- /dev/null
+++ b/PARAM/cation_param.parm
@@ -0,0 +1,11 @@
+ 23.0  1.868 
+ 24.0  0.793
+ 39.0  2.650
+ 40.0  1.713
+ 35.5  2.470
+7
+-1381.06 282.2  -203.7  3.0   0.1      3.5   0.12
+-1426.87 -98.92 -2454.0 3.0   0.1      3.5   0.02
+0        -70.79 1499.0  3.0   0.1      3.5   0.05
+0        -93.13 1975.0  3.0   0.1      3.5   0.01
+wc        wdip wquad1 wquad2 epsilonLJ r0lJ sigm

diff --git a/source/unres/MD.f90 b/source/unres/MD.f90
index 82ecbe8..1af6b18 100644 (file)
--- a/source/unres/MD.f90
+++ b/source/unres/MD.f90
@@ -652,7 +652,7 @@
 !      include 'COMMON.INTERACT'
 !      include 'COMMON.MD'
 !      include 'COMMON.IOUNITS'
-      real(kind=8) :: KE_total
+      real(kind=8) :: KE_total,mscab
                                                              
       integer :: i,j,k,iti,mnum
       real(kind=8) :: KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3),&
@@ -673,6 +673,7 @@
       enddo
       do i=nnt,nct-1
        mnum=molnum(i)
+       if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
 !        write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3) 
         do j=1,3
           v(j)=incr(j)+0.5d0*d_t(j,i)
@@ -693,6 +694,11 @@
       do i=nnt,nct
          mnum=molnum(i)
         iti=iabs(itype(i,mnum))
+        if (mnum.eq.5) then
+         mscab=0.0d0
+        else
+         mscab=msc(iti,mnum)
+        endif
 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
          if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
           .or.(mnum.eq.5)) then
@@ -706,7 +712,7 @@
         endif
 !        write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3) 
 !        write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3) 
-        KEt_sc=KEt_sc+msc(iti,mnum)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))            
+        KEt_sc=KEt_sc+mscab*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))            
         vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3)
         do j=1,3
           incr(j)=incr(j)+d_t(j,i)
@@ -977,7 +983,8 @@
           enddo
         enddo
         do i=nnt,nct
-          if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
+          mnum=molnum(i)
+          if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.ne.5) then
             do j=1,3
               ind=ind+1
               v_work(ind)=d_t(j,i+nres)
@@ -2209,7 +2216,8 @@
           enddo
         endif
 ! Side chains
-        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
+        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.&
+        molnum(i).ne.5) then
           do j=1,3 
             epdriftij= &
              dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
@@ -2911,6 +2919,7 @@
       do i=nnt,nct
 !        if (itype(i,1).eq.10) then
          mnum=molnum(i)
+          if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
         iti=iabs(itype(i,mnum))
 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
          if (itype(i,1).eq.10 .or. itype(i,mnum).eq.ntyp1_molec(mnum)&
@@ -3742,6 +3751,7 @@
         M_PEP=0.0d0
         do i=nnt,nct-1
           mnum=molnum(i)
+          if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
           if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
           M_PEP=M_PEP+mp(mnum)
           do j=1,3
@@ -3754,6 +3764,7 @@
         M_SC=0.0d0                             
         do i=nnt,nct
            mnum=molnum(i)
+           if (mnum.eq.5) cycle
            iti=iabs(itype(i,mnum))              
           M_SC=M_SC+msc(iabs(iti),mnum)
            inres=i+nres
@@ -3767,6 +3778,7 @@
        
         do i=nnt,nct-1
            mnum=molnum(i)
+          if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
           do j=1,3
             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
           enddo
@@ -3781,6 +3793,7 @@
        do i=nnt,nct    
            mnum=molnum(i)
            iti=iabs(itype(i,mnum))
+           if (mnum.eq.5) cycle
            inres=i+nres
            do j=1,3
              pr(j)=c(j,inres)-cm(j)        
@@ -3931,7 +3944,7 @@
 !       include 'COMMON.INTERACT'
 !       include 'COMMON.IOUNITS'
 !       include 'COMMON.NAMES'
-      
+      real(kind=8) :: mscab
       real(kind=8),dimension(3) :: L,cm,pr,vp,vrot,incr,v,pp
       integer :: iti,inres,i,j,mnum
 !  Calculate the angular momentum
@@ -3943,6 +3956,7 @@
        enddo                  
        do i=nnt,nct-1
           mnum=molnum(i)
+          if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
           do j=1,3
             pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j)
           enddo
@@ -3972,6 +3986,11 @@
           mnum=molnum(i)
          iti=iabs(itype(i,mnum))
          inres=i+nres
+        if (mnum.eq.5) then
+         mscab=0.0d0
+        else
+         mscab=msc(iti,mnum)
+        endif
          do j=1,3
            pr(j)=c(j,inres)-cm(j)          
          enddo
@@ -3989,7 +4008,7 @@
 !         write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3),
 !     &     " v",(v(j),j=1,3)," vp",(vp(j),j=1,3)
          do j=1,3
-            L(j)=L(j)+msc(iabs(iti),mnum)*vp(j)
+            L(j)=L(j)+mscab*vp(j)
          enddo
 !         write (iout,*) "L",(l(j),j=1,3)
          if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
@@ -4033,13 +4052,18 @@
        summas=0.0d0
        do i=nnt,nct
          mnum=molnum(i)
+         if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
          if (i.lt.nct) then
            summas=summas+mp(mnum)
            do j=1,3
              vcm(j)=vcm(j)+mp(mnum)*(vv(j)+0.5d0*d_t(j,i))
            enddo
          endif
+         if (mnum.ne.5) then 
          amas=msc(iabs(itype(i,mnum)),mnum)
+         else
+         amas=0.0d0
+         endif
          summas=summas+amas                     
          if (itype(i,mnum).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
           .and.(mnum.ne.5)) then

diff --git a/source/unres/REMD.f90 b/source/unres/REMD.f90
index 1e294ad..c85328c 100644 (file)
--- a/source/unres/REMD.f90
+++ b/source/unres/REMD.f90
@@ -47,7 +47,7 @@
 #ifdef CHECK5DSOL
        real(kind=8) :: rscheck(dimen),rsold(dimen)
 #endif
-       logical :: lprn = .true.
+       logical :: lprn = .false.
 !el       common /cipiszcze/ itime
        itime = itt_comm
 
@@ -62,14 +62,16 @@
         enddo
         write (iout,*) "Potential forces sidechain"
         do i=nnt,nct
-          if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) &
+          mnum=molnum(i)
+          if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1.and.mnum.eq.5)&
             write (iout,'(i5,3e15.5,5x,3e15.5)') i,(-gxcart(j,i),j=1,3)
         enddo
       endif
        do j=1,3
          ind=1
          do i=nnt,nct
-           if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1)then
+           mnum=molnum(i)     
+           if (itype(i,1).eq.10 .or. itype(i,1).eq.ntyp1 .or. mnum.eq.5)then
            rs(ind)=-gcart(j,i)-gxcart(j,i)
             ind=ind+1
           else
@@ -192,7 +194,8 @@
       if (lprn) write (iout,*) "Potential forces sidechain"
       do i=nnt,nct
           mnum=molnum(i)
-        if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
+        if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+          .and.mnum.ne.5) then
           if (lprn) write (iout,'(i5,3e15.5,5x,3e15.5)') &
              i,(-gxcart(j,i),j=1,3)
           do j=1,3
@@ -217,7 +220,8 @@
       do i=nnt,nct
        mnum=molnum(i)
 !        if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
-         if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
+         if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)&
+          .and.mnum.ne.5) then
           do j=1,3
             ind=ind+1
             d_a(j,i+nres)=d_a_work(ind)
@@ -311,7 +315,7 @@
       real(kind=8),dimension(8*6*nres) :: work !(8*maxres6)
       integer,dimension(6*nres) :: iwork       !(maxres6) maxres6=6*maxres
 !el      common /przechowalnia/ Gcopy,Ghalf
-      real(kind=8) :: coeff
+      real(kind=8) :: coeff,mscab
       integer :: i,j,ind,ind1,k,l,ii,jj,m,m1,ii1,iti,nres2,ierr,nind,mark
       real(kind=8) :: ip4
       integer :: iz,mnum
@@ -374,7 +378,9 @@
        
         ind=1
         do i=nnt,nct
-       if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1) then
+        mnum=molnum(i)
+       if (iabs(itype(i,1)).ne.10 .and.iabs(itype((i))).ne.ntyp1 &
+          .and.mnum.eq.5) then
           DU1(ind)=-isc(iabs(itype(i,1)))
           DU1(ind+1)=0.0d0
          ind=ind+2
@@ -513,6 +519,8 @@
         ind=ind+1
         ind1=ind1+1
         coeff=0.25d0*IP(mnum)
+        if (mnum.eq.5) mp(mnum)=msc(itype(i,mnum),mnum)
+        print *,i,coeff,mp(mnum)
         massvec(ind1)=mp(mnum)
         Gmat(ind,ind)=coeff
         print *,"i",mp(mnum)
@@ -539,7 +547,13 @@
         ii = ind+m
         mnum=molnum(i)
         iti=itype(i,mnum)
-        massvec(ii)=msc(iabs(iti),mnum)
+        if (mnum.eq.5) then
+        mscab=0.0
+        else
+        mscab=msc(iabs(iti),mnum)
+        endif
+        massvec(ii)=mscab
+
         if (iti.ne.10 .and. iti.ne.ntyp1_molec(mnum) .and. mnum.ne.5) then
           ind1=ind1+1
           ii1= ind1+m1

diff --git a/source/unres/compare.F90 b/source/unres/compare.F90
index 8e6c0af..50972e1 100644 (file)
--- a/source/unres/compare.F90
+++ b/source/unres/compare.F90
@@ -818,7 +818,7 @@
 !     &             obr,non_conv)
 !        rms=dsqrt(rms)
         call rmsd(rms)
-        print *,"before contact"
+!        print *,"before contact"
 !elte(iout,*) "rms_nacc before contact"
         call contact(.false.,ncont,icont,co)
         frac=contact_fract(ncont,ncont_ref,icont,icont_ref)

diff --git a/source/unres/control.F90 b/source/unres/control.F90
index 7ecc757..fbd2dfb 100644 (file)
--- a/source/unres/control.F90
+++ b/source/unres/control.F90
@@ -242,6 +242,8 @@
 
       iliptranpar=60
       itube=61
+!     IONS
+      iion=401
 #if defined(WHAM_RUN) || defined(CLUSTER)
 !
 ! setting the mpi variables for WHAM
@@ -1921,7 +1923,7 @@
 #ifdef WHAM_RUN
         if (itype(i,1).ne.10) then
 #else
-        if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum)) then
+        if (itype(i,1).ne.10 .and. itype(i,mnum).ne.ntyp1_molec(mnum) .and. mnum.ne.5) then
 #endif
          nside=nside+1
           ialph(i,1)=nvar+nside

diff --git a/source/unres/data/energy_data.f90 b/source/unres/data/energy_data.f90
index baf7f55..8807974 100644 (file)
--- a/source/unres/data/energy_data.f90
+++ b/source/unres/data/energy_data.f90
@@ -66,7 +66,8 @@
        wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,&
        wturn6,wvdwpp,wliptran,wshield,lipscale,wtube, &
        wbond_nucl,wang_nucl,wcorr_nucl,wcorr3_nucl,welpp,wtor_nucl,&
-       wtor_d_nucl,welsb,wsbloc,wvdwsb,welpsb,wvdwpp_nucl,wvdwpsb
+       wtor_d_nucl,welsb,wsbloc,wvdwsb,welpsb,wvdwpp_nucl,wvdwpsb,wcatprot,&
+       wcatcat
 #ifdef CLUSTER
       real(kind=8) :: scalscp
 #endif
@@ -336,7 +337,7 @@
 ! AFM
        real(kind=8) :: distafminit,forceAFMconst,velAFMconst
       integer :: afmend,afmbeg
-
+      real(kind=8),dimension(:,:), allocatable :: catprm
 
 
 

diff --git a/source/unres/data/io_units.f90 b/source/unres/data/io_units.f90
index 5e9e9b7..439ed89 100644 (file)
--- a/source/unres/data/io_units.f90
+++ b/source/unres/data/io_units.f90
@@ -16,7 +16,8 @@
        ientout,izs1,isecpred,ibond,irest2,iifrag,icart,irest1,isccor,&
        ithep_pdb,irotam_pdb,iliptranpar,itube,   &
        ibond_nucl,ithep_nucl,irotam_nucl,itorp_nucl,itordp_nucl,     &
-       isidep_nucl,iscpp_nucl      
+       isidep_nucl,iscpp_nucl,   &
+       iion      
 #ifdef WHAM_RUN
 ! el wham iounits
       integer :: isidep1,ihist,iweight,izsc,idistr
@@ -47,7 +48,8 @@
        fouriername,elename,sidename,scpname,sccorname,patname,&
        thetname_pdb,rotname_pdb,liptranname,tubename,         &
        bondname_nucl,thetname_nucl,rotname_nucl,torname_nucl, &
-       tordname_nucl,sidename_nucl,scpname_nucl
+       tordname_nucl,sidename_nucl,scpname_nucl,              &
+       ionname
 !-----------------------------------------------------------------------
 ! INP    - main input file
 ! IOUT   - list file

diff --git a/source/unres/data/names.f90 b/source/unres/data/names.f90
index cb3f2f9..1d82366 100644 (file)
--- a/source/unres/data/names.f90
+++ b/source/unres/data/names.f90
@@ -59,7 +59,7 @@
 !-----------------------------------------------------------------------------
 !-----------------------------------------------------------------------------
 ! Number of energy components
-      integer,parameter :: n_ene=38
+      integer,parameter :: n_ene=42
       integer :: n_ene2=2*n_ene
 !-----------------------------------------------------------------------------
 ! common.names
@@ -86,7 +86,8 @@
         "ELT       ","          ","          ","ETUBE     ",&
         "EVDWPP    ","EESPP     ","EVDWPSB   ","EESPSB    ","EVDWSB    ",&
         "EESSB     ","ESTR      ","EBE       ","ESBLOC    ","ETORS     ",&
-        "ETORSD    ","ECORR     ","ECORR3    " /)
+        "ETORSD    ","ECORR     ","ECORR3    ","NULL      ","NULL      ",&
+        "ECATPROT  ","ECATCAT   "/)
 
       character(len=10),dimension(n_ene) :: wname = &
       (/"WSC       ","WSCP      ","WELEC"    ,"WCORR      ","WCORR5    ","WCORR6    ","WEL_LOC   ",&
@@ -95,12 +96,13 @@
         "WLT       ","          ","         ","WTUBE      " ,&
         "WVDWPP    ","WELPP     ","WVDWPSB  ","WELPSB     ","WVDWSB    ",&
         "WELSB     ","WBOND     ","WANG     ","WSBLOC     ","WTOR      ",&
-        "WTORD     ","WCORR     ","WCORR3   "/)
+        "WTORD     ","WCORR     ","WCORR3   ","WNULL      ","WNULL     ",&
+        "WCATPROT  ","WCATCAT   "/)
 
       integer :: nprint_ene = 21
       integer,dimension(n_ene) :: print_order = &
          (/1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,16,15,17,20,21,22,23,24,25,&
-         26,27,28,29,30,31,32,33,34,35,36,37,38/)
+         26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42/)
 
       character(len=1), dimension(2) :: sugartyp = (/'D',' '/)
 !#endif

diff --git a/source/unres/energy.f90 b/source/unres/energy.f90
index 964624f..c92f4bc 100644 (file)
--- a/source/unres/energy.f90
+++ b/source/unres/energy.f90
@@ -131,7 +131,8 @@
         gsblocx,gradcorr_nucl,gradxorr_nucl,gradcorr3_nucl,gradxorr3_nucl,&
         gvdwpp_nucl
 !------------------------------IONS GRADIENT
-        real(kind=8),dimension(:,:),allocatable  ::  gradcatcat
+        real(kind=8),dimension(:,:),allocatable  ::  gradcatcat, &
+          gradpepcat,gradpepcatx
 !      real(kind=8),dimension(:,:),allocatable :: gloc,gloc_x !(maxvar,2)
       real(kind=8),dimension(:,:),allocatable :: gel_loc,gel_loc_long,&
         gcorr3_turn,gcorr4_turn,gcorr6_turn,gradb,gradbx !(3,maxres)
@@ -238,6 +239,9 @@
       real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,&
                       ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
                       ecorr3_nucl
+! energies for ions 
+      real(kind=8) :: ecation_prot,ecationcation
+
 #ifdef MPI      
       real(kind=8) :: weights_(n_ene) !,time_Bcast,time_Bcastw
 ! shielding effect varibles for MPI
@@ -290,6 +294,10 @@
           weights_(36)=wtor_d_nucl
           weights_(37)=wcorr_nucl
           weights_(38)=wcorr3_nucl
+          weights(41)=wcatcat
+          weights(42)=wcatprot
+!          wcatcat= weights(41)
+!          wcatprot=weights(42)
 
 ! FG Master broadcasts the WEIGHTS_ array
           call MPI_Bcast(weights_(1),n_ene,&
@@ -330,7 +338,8 @@
           wtor_d_nucl =weights(36)
           wcorr_nucl  =weights(37)
           wcorr3_nucl =weights(38)
-
+          wcatcat= weights(41)
+          wcatprot=weights(42)
         endif
         time_Bcast=time_Bcast+MPI_Wtime()-time00
         time_Bcastw=time_Bcastw+MPI_Wtime()-time00
@@ -590,7 +599,9 @@
       call epsb(evdwpsb,eelpsb)
       call esb(esbloc)
       call multibody_hb_nucl(ecorr_nucl,ecorr3_nucl,n_corr,n_corr1)
-
+      call ecatcat(ecationcation)
+      call ecat_prot(ecation_prot)
+!      call ecatcat(ecationcation)
 !      print *,"after ebend", ebe_nucl
 #ifdef TIMING
       time_enecalc=time_enecalc+MPI_Wtime()-time00
@@ -655,6 +666,8 @@
 !    Here are the energies showed per procesor if the are more processors 
 !    per molecule then we sum it up in sum_energy subroutine 
 !      print *," Processor",myrank," calls SUM_ENERGY"
+      energia(41)=ecation_prot
+      energia(42)=ecationcation
       call sum_energy(energia,.true.)
       if (dyn_ss) call dyn_set_nss
 !      print *," Processor",myrank," left SUM_ENERGY"
@@ -697,7 +710,7 @@
       real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,&
                       ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
                       ecorr3_nucl
-
+      real(kind=8) :: ecation_prot,ecationcation
       integer :: i
 #ifdef MPI
       integer :: ierr
@@ -774,6 +787,10 @@
       etors_d_nucl=energia(36)
       ecorr_nucl=energia(37)
       ecorr3_nucl=energia(38)
+      ecation_prot=energia(41)
+      ecationcation=energia(42)
+!      energia(41)=ecation_prot
+!      energia(42)=ecationcation
 
 
 #ifdef SPLITELE
@@ -787,7 +804,8 @@
        +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
        +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
        +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
-       +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl
+       +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
+       +wcatprot*ecation_prot+wcatcat*ecationcation
 #else
       etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) &
        +wang*ebe+wtor*etors+wscloc*escloc &
@@ -799,7 +817,8 @@
        +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
        +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
        +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
-       +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl
+       +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
+       +wcatprot*ecation_prot+wcatcat*ecationcation
 #endif
       energia(0)=etot
 ! detecting NaNQ
@@ -927,6 +946,7 @@
       real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,&
                       ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
                       ecorr3_nucl
+      real(kind=8) :: ecation_prot,ecationcation
 
       etot=energia(0)
       evdw=energia(1)
@@ -973,6 +993,8 @@
       etors_d_nucl=energia(36)
       ecorr_nucl=energia(37)
       ecorr3_nucl=energia(38)
+      ecation_prot=energia(41)
+      ecationcation=energia(42)
 
 #ifdef SPLITELE
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,&
@@ -987,7 +1009,7 @@
         evdwpp,wvdwpp_nucl,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb,&
         evdwsb,wvdwsb,eelsb,welsb,esbloc,wsbloc,etors_nucl,wtor_nucl,&
         etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,&
-        ecorr3_nucl,wcorr3_nucl, &
+        ecorr3_nucl,wcorr3_nucl,ecation_prot,wcatprot,ecationcation,wcatcat, &
         etot
    10 format (/'Virtual-chain energies:'// &
        'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ &
@@ -1029,6 +1051,8 @@
        'ETORSD_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(double torsional)'/ &
        'ECORR_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(multibody 4th order)'/ &
        'ECORR3_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(multibody 3th order)'/ &
+       'ECATPROT=',1pE16.6,' WEIGHT=',1pD16.6,'(cation prot)'/ &
+       'ECATCAT=',1pE16.6,' WEIGHT=',1pD16.6,'(cation cation)'/ &
        'ETOT=  ',1pE16.6,' (total)')
 #else
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,&
@@ -1043,7 +1067,7 @@
         evdwpp,wvdwpp_nucl,eespp,welpp,evdwpsb,wvdwpsb,eelpsb,welpsb&
         evdwsb,wvdwsb,eelsb,welsb,esbloc,wsbloc,etors_nucl,wtor_nucl&
         etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,&
-        ecorr3_nucl,wcorr3_nucl, &
+        ecorr3_nucl,wcorr3_nucl,ecation_prot,wcatprot,ecationcation,wcatcat,  &
         etot
    10 format (/'Virtual-chain energies:'// &
        'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ &
@@ -1084,6 +1108,8 @@
        'ETORSD_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(double torsional)'/ &
        'ECORR_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(multibody 4th order)'/ &
        'ECORR3_nucl=',1pE16.6,' WEIGHT=',1pD16.6,'(multibody 3th order)'/ &
+       'ECATPROT=',1pE16.6,' WEIGHT=',1pD16.6,'(cation prot)'/ &
+       'ECATCAT=',1pE16.6,' WEIGHT=',1pD16.6,'(cation cation)'/ &
        'ETOT=  ',1pE16.6,' (total)')
 #endif
       return
@@ -10448,7 +10474,10 @@
                      wvdwpsb*(gvdwpsb(j,i)+gvdwpsb1(j,i))+ &
                      wvdwsb*gvdwsbc(j,i)+welsb*gelsbc(j,i)+ &
                      wcorr_nucl*gradcorr_nucl(j,i)&
-                     +wcorr3_nucl*gradcorr3_nucl(j,i)
+                     +wcorr3_nucl*gradcorr3_nucl(j,i)+&
+                     wcatprot* gradpepcat(j,i)+ &
+                     wcatcat*gradcatcat(j,i)
+
 
         enddo
       enddo 
@@ -10475,8 +10504,10 @@
                      +wvdwpp_nucl*gvdwpp_nucl(j,i)+welpp*gelpp(j,i)+ &
                      wvdwpsb*(gvdwpsb(j,i)+gvdwpsb1(j,i))+ &
                      wvdwsb*gvdwsbc(j,i)+welsb*gelsbc(j,i)+ &
-                     wcorr_nucl*gradcorr_nucl(j,i)
-                     +wcorr3_nucl*gradcorr3_nucl(j,i)
+                     wcorr_nucl*gradcorr_nucl(j,i) &
+                     +wcorr3_nucl*gradcorr3_nucl(j,i) +&
+                     wcatprot* gradpepcat(j,i)+ &
+                     wcatcat*gradcatcat(j,i)
         enddo
       enddo 
 #endif
@@ -10734,7 +10765,8 @@
                        +welsb*gelsbx(j,i) &
                        +wcorr_nucl*gradxorr_nucl(j,i)&
                        +wcorr3_nucl*gradxorr3_nucl(j,i) &
-                       +wsbloc*gsblocx(j,i)
+                       +wsbloc*gsblocx(j,i) &
+                       +wcatprot* gradpepcatx(j,i)
         enddo
       enddo 
 #ifdef DEBUG
@@ -16169,9 +16201,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
 #ifdef TIMING
       time01=MPI_Wtime()
 #endif
-       print *,"gcart_two",gxcart(2,2),gradx(2,2,icg)
+!       print *,"gcart_two",gxcart(2,2),gradx(2,2,icg)
       call int_to_cart
-             print *,"gcart_two",gxcart(2,2),gradx(2,2,icg)
+!             print *,"gcart_two",gxcart(2,2),gradx(2,2,icg)
 
 #ifdef TIMING
             time_inttocart=time_inttocart+MPI_Wtime()-time01
@@ -16338,6 +16370,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
               gelsbx(j,i)=0.0d0
               gsbloc(j,i)=0.0d0
               gsblocx(j,i)=0.0d0
+              gradpepcat(j,i)=0.0d0
+              gradpepcatx(j,i)=0.0d0
+              gradcatcat(j,i)=0.0d0
             enddo
              enddo
             do i=0,nres
@@ -19717,7 +19752,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
       allocate(gradcorr3_nucl(3,-1:nres))
       allocate(gradxorr3_nucl(3,-1:nres))
       allocate(gvdwpp_nucl(3,-1:nres))
-
+      allocate(gradpepcat(3,-1:nres))
+      allocate(gradpepcatx(3,-1:nres))
+      allocate(gradcatcat(3,-1:nres))
 !(3,maxres)
       allocate(grad_shield_side(3,50,nres))
       allocate(grad_shield_loc(3,50,nres))
@@ -21622,12 +21659,13 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         epscalc=0.05
         r06 = rcat0**6
         r012 = r06**2
-        k0 = 332*(2*2)/80
+        k0 = 332.0*(2.0*2.0)/80.0
         itmp=0
         do i=1,4
         itmp=itmp+nres_molec(i)
         enddo
         do i=itmp+1,itmp+nres_molec(i)-1
+       
         xi=c(1,i)
         yi=c(2,i)
         zi=c(3,i)
@@ -21638,7 +21676,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
           zi=mod(zi,boxzsize)
           if (zi.lt.0) zi=zi+boxzsize
 
-          do j=i+1,itmp+nres_molec(i)
+          do j=i+1,itmp+nres_molec(5)
+!           print *,i,j,'catcat'
            xj=c(1,j)
            yj=c(2,j)
            zj=c(3,j)
@@ -21701,8 +21740,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
         enddo
         do k=1,3
           gg(k) = dEvan1Cmcat(k)+dEvan2Cmcat(k)+dEeleccat(k)
-          gradcatcat(k,i)=gradcatcat(k,i)+gg(k)
-          gradcatcat(k,j)=gradcatcat(k,j)-gg(k)
+          gradcatcat(k,i)=gradcatcat(k,i)-gg(k)
+          gradcatcat(k,j)=gradcatcat(k,j)+gg(k)
         enddo
 
         ecationcation=ecationcation+Evan1cat+Evan2cat+Eeleccat
@@ -21712,20 +21751,50 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
        end subroutine ecatcat
 !---------------------------------------------------------------------------
        subroutine ecat_prot(ecation_prot)
-       integer i,j,k,subchap,itmp
+       integer i,j,k,subchap,itmp,inum
         real(kind=8) :: xi,yi,zi,xj,yj,zj,ract,rcat0,epscalc,r06,r012,&
         r7,r4,ecationcation
         real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
-        dist_init,dist_temp,ecation_prot
+        dist_init,dist_temp,ecation_prot,rcal,rocal,   &
+        Evan1,Evan2,EC,cm1mag,DASGL,delta,r0p,Epepcat, &
+        catl,cml,calpl, Etotal_p, Etotal_m,rtab,wdip,wmodquad,wquad1, &
+        wquad2,wvan1,E1,E2,wconst,wvan2,rcpm,dcmag,sin2thet,sinthet,  &
+        costhet,v1m,v2m,wh2o,wc,rsecp,Ir,Irsecp,Irthrp,Irfourp,Irfiftp,&
+        Irsistp,Irseven,Irtwelv,Irthir,dE1dr,dE2dr,dEdcos,wquad2p,opt, &
+        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
         real(kind=8),dimension(3) ::dEvan1Cmcat,dEvan2Cmcat,dEeleccat,&
-        gg,r
+        gg,r,EtotalCat,dEtotalCm,dEtotalCalp,dEvan1Cm,dEvan2Cm, &
+        dEtotalpep,dEtotalcat_num,dEddci,dEtotalcm_num,dEtotalcalp_num, &
+        tab1,tab2,tab3,diff,cm1,sc,p,tcat,talp,cm,drcp,drcp_norm,vcat,  &
+        v1,v2,v3,myd_norm,dx,vcm,valpha,drdpep,dcosdpep,dcosddci,dEdpep,&
+        dEcCat,dEdipCm,dEdipCalp,dEquad1Cat,dEquad1Cm,dEquad1Calp,      &
+        dEquad2Cat,dEquad2Cm,dEquad2Calpd,Evan1Cat,dEvan1Calp,dEvan2Cat,&
+        dEvan2Calp,dEtotalCat,dscvec,dEcCm,dEcCalp,dEdipCat,dEquad2Calp,&
+        dEvan1Cat
+        real(kind=8),dimension(6) :: vcatprm
+        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)
         enddo
         do i=1,nres_molec(1)-1  ! loop over all peptide groups needs parralelization
+!         cycle
          if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle ! leave dummy atoms
         xi=0.5d0*(c(1,i)+c(1,i+1))
         yi=0.5d0*(c(2,i)+c(2,i+1))
@@ -21778,10 +21847,446 @@ write(iout,*) 'Calling CHECK_ECARTIN else.'
           yj=yj_safe-yi
           zj=zj_safe-zi
        endif
+!       enddo
+!       enddo
+       rcpm = sqrt(xj**2+yj**2+zj**2)
+       drcp_norm(1)=xj/rcpm
+       drcp_norm(2)=yj/rcpm
+       drcp_norm(3)=zj/rcpm
+       dcmag=0.0
+       do k=1,3
+       dcmag=dcmag+dc(k,i)**2
        enddo
+       dcmag=dsqrt(dcmag)
+       do k=1,3
+         myd_norm(k)=dc(k,i)/dcmag
        enddo
+        costhet=drcp_norm(1)*myd_norm(1)+drcp_norm(2)*myd_norm(2)+&
+        drcp_norm(3)*myd_norm(3)
+        rsecp = rcpm**2
+        Ir = 1.0d0/rcpm
+        Irsecp = 1.0d0/rsecp
+        Irthrp = Irsecp/rcpm
+        Irfourp = Irthrp/rcpm
+        Irfiftp = Irfourp/rcpm
+        Irsistp=Irfiftp/rcpm
+        Irseven=Irsistp/rcpm
+        Irtwelv=Irsistp*Irsistp
+        Irthir=Irtwelv/rcpm
+        sin2thet = (1-costhet*costhet)
+        sinthet=sqrt(sin2thet)
+        E1 = wdip*Irsecp*costhet+(wmodquad*Irfourp+wquad1*Irthrp)&
+             *sin2thet
+        E2 = -wquad1*Irthrp*wquad2+wvan1*(wvan2**12*Irtwelv-&
+             2*wvan2**6*Irsistp)
+        ecation_prot = ecation_prot+E1+E2
+        dE1dr = -2*costhet*wdip*Irthrp-& 
+         (4*wmodquad*Irfiftp+3*wquad1*Irfourp)*sin2thet
+        dE2dr = 3*wquad1*wquad2*Irfourp-     &
+          12*wvan1*wvan2**6*(wvan2**6*Irthir-Irseven)
+        dEdcos = wdip*Irsecp-2*(wmodquad*Irfourp+wquad1*Irthrp)*costhet
+        do k=1,3
+          drdpep(k) = -drcp_norm(k)
+          dcosdpep(k) = Ir*(costhet*drcp_norm(k)-myd_norm(k))
+          dcosddci(k) = drcp_norm(k)/dcmag-costhet*myd_norm(k)/dcmag
+          dEdpep(k) = (dE1dr+dE2dr)*drdpep(k)+dEdcos*dcosdpep(k)
+          dEddci(k) = dEdcos*dcosddci(k)
+        enddo
+        do k=1,3
+        gradpepcat(k,i)=gradpepcat(k,i)+0.5D0*dEdpep(k)-dEddci(k)
+        gradpepcat(k,i+1)=gradpepcat(k,i+1)+0.5D0*dEdpep(k)+dEddci(k)
+        gradpepcat(k,j)=gradpepcat(k,j)-dEdpep(k)
+        enddo
+       enddo ! j
+       enddo ! i
+!------------------------------------------sidechains
+        do i=1,nres_molec(1)
+         if ((itype(i,1).eq.ntyp1)) cycle ! leave dummy atoms
+!         cycle
+!        print *,i,ecation_prot
+        xi=(c(1,i+nres))
+        yi=(c(2,i+nres))
+        zi=(c(3,i+nres))
+          xi=mod(xi,boxxsize)
+          if (xi.lt.0) xi=xi+boxxsize
+          yi=mod(yi,boxysize)
+          if (yi.lt.0) yi=yi+boxysize
+          zi=mod(zi,boxzsize)
+          if (zi.lt.0) zi=zi+boxzsize
+          do k=1,3
+            cm1(k)=dc(k,i+nres)
+          enddo
+           cm1mag=sqrt(cm1(1)**2+cm1(2)**2+cm1(3)**2)
+         do j=itmp+1,itmp+nres_molec(5)
+           xj=c(1,j)
+           yj=c(2,j)
+           zj=c(3,j)
+          xj=dmod(xj,boxxsize)
+          if (xj.lt.0) xj=xj+boxxsize
+          yj=dmod(yj,boxysize)
+          if (yj.lt.0) yj=yj+boxysize
+          zj=dmod(zj,boxzsize)
+          if (zj.lt.0) zj=zj+boxzsize
+      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+      xj_safe=xj
+      yj_safe=yj
+      zj_safe=zj
+      subchap=0
+      do xshift=-1,1
+      do yshift=-1,1
+      do zshift=-1,1
+          xj=xj_safe+xshift*boxxsize
+          yj=yj_safe+yshift*boxysize
+          zj=zj_safe+zshift*boxzsize
+          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+          if(dist_temp.lt.dist_init) then
+            dist_init=dist_temp
+            xj_temp=xj
+            yj_temp=yj
+            zj_temp=zj
+            subchap=1
+          endif
+       enddo
+       enddo
+       enddo
+       if (subchap.eq.1) then
+          xj=xj_temp-xi
+          yj=yj_temp-yi
+          zj=zj_temp-zi
+       else
+          xj=xj_safe-xi
+          yj=yj_safe-yi
+          zj=zj_safe-zi
+       endif
+!       enddo
+!       enddo
+         if(itype(i,1).eq.15.or.itype(i,1).eq.16) then
+            if(itype(i,1).eq.16) then
+            inum=1
+            else
+            inum=2
+            endif
+            do k=1,6
+            vcatprm(k)=catprm(k,inum)
+            enddo
+            dASGL=catprm(7,inum)
+             do k=1,3
+                vcm(k)=(cm1(k)/cm1mag)*dASGL+c(k,i+nres)
+                valpha(k)=c(k,i)
+                vcat(k)=c(k,j)
+              enddo
+                      do k=1,3
+          dx(k) = vcat(k)-vcm(k)
+        enddo
+        do k=1,3
+          v1(k)=(vcm(k)-valpha(k))
+          v2(k)=(vcat(k)-valpha(k))
+        enddo
+        v1m = sqrt(v1(1)**2+v1(2)**2+v1(3)**2)
+        v2m = sqrt(v2(1)**2+v2(2)**2+v2(3)**2)
+        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 calcium with ASP/GLU
+        wh2o=78
+        wc = vcatprm(1)
+        wc=wc/wh2o
+        wdip =vcatprm(2)
+        wdip=wdip/wh2o
+        wquad1 =vcatprm(3)
+        wquad1=wquad1/wh2o
+        wquad2 = vcatprm(4)
+        wquad2=wquad2/wh2o
+        wquad2p = 1-wquad2
+        wvan1 = vcatprm(5)
+        wvan2 =vcatprm(6)
+        opt = dx(1)**2+dx(2)**2
+        rsecp = opt+dx(3)**2
+        rs = sqrt(rsecp)
+        rthrp = rsecp*rs
+        rfourp = rthrp*rs
+        rsixp = rfourp*rsecp
+        reight=rsixp*rsecp
+        Ir = 1.0d0/rs
+        Irsecp = 1/rsecp
+        Irthrp = Irsecp/rs
+        Irfourp = Irthrp/rs
+        Irsixp = 1/rsixp
+        Ireight=1/reight
+        Irtw=Irsixp*Irsixp
+        Irthir=Irtw/rs
+        Irfourt=Irthir/rs
+        opt1 = (4*rs*dx(3)*wdip)
+        opt2 = 6*rsecp*wquad1*opt
+        opt3 = wquad1*wquad2p*Irsixp
+        opt4 = (wvan1*wvan2**12)
+        opt5 = opt4*12*Irfourt
+        opt6 = 2*wvan1*wvan2**6
+        opt7 = 6*opt6*Ireight
+        opt8 = wdip/v1m
+        opt10 = wdip/v2m
+        opt11 = (rsecp*v2m)**2
+        opt12 = (rsecp*v1m)**2
+        opt14 = (v1m*v2m*rsecp)**2
+        opt15 = -wquad1/v2m**2
+        opt16 = (rthrp*(v1m*v2m)**2)**2
+        opt17 = (v1m**2*rthrp)**2
+        opt18 = -wquad1/rthrp
+        opt19 = (v1m**2*v2m**2)**2
+        Ec = wc*Ir
+        do k=1,3
+          dEcCat(k) = -(dx(k)*wc)*Irthrp
+          dEcCm(k)=(dx(k)*wc)*Irthrp
+          dEcCalp(k)=0.0d0
+        enddo
+        Edip=opt8*(v1dpv2)/(rsecp*v2m)
+        do k=1,3
+          dEdipCat(k)=opt8*(v1(k)*rsecp*v2m-((v2(k)/v2m &
+                     *rsecp+2*dx(k)*v2m)*v1dpv2))/opt11
+          dEdipCm(k)=opt10*(v2(k)*rsecp*v1m-((v1(k)/v1m &
+                    *rsecp-2*dx(k)*v1m)*v1dpv2))/opt12
+          dEdipCalp(k)=wdip*((-v1(k)-v2(k))*rsecp*v1m &
+                      *v2m-(-v1(k)/v1m*v2m*rsecp-v2(k)/v2m*v1m*rsecp) &
+                      *v1dpv2)/opt14
+        enddo
+        Equad1=-wquad1*v1dpv2**2/(rthrp*(v1m*v2m)**2)
+        do k=1,3
+          dEquad1Cat(k)=-wquad1*(2*v1(k)*v1dpv2*(rthrp* &
+                       (v1m*v2m)**2)-(3*dx(k)*rs*(v1m*v2m)**2+2*v1m*2* &
+                       v2(k)*1/2*1/v2m*v1m*v2m*rthrp)*v1dpv2**2)/opt16
+          dEquad1Cm(k)=-wquad1*(2*v2(k)*v1dpv2*(rthrp* &
+                      (v1m*v2m)**2)-(-3*dx(k)*rs*(v1m*v2m)**2+2*v2m*2* &
+                      v1(k)*1/2*1/v1m*v2m*v1m*rthrp)*v1dpv2**2)/opt16
+          dEquad1Calp(k)=opt18*(2*(-v1(k)-v2(k))*v1dpv2* &
+                        v1m**2*v2m**2-(-2*v1(k)*v2m**2-2*v2(k)*v1m**2)* &
+                        v1dpv2**2)/opt19
+        enddo
+        Equad2=wquad1*wquad2p*Irthrp
+        do k=1,3
+          dEquad2Cat(k)=-3*dx(k)*rs*opt3
+          dEquad2Cm(k)=3*dx(k)*rs*opt3
+          dEquad2Calp(k)=0.0d0
+        enddo
+        Evan1=opt4*Irtw
+        do k=1,3
+          dEvan1Cat(k)=-dx(k)*opt5
+          dEvan1Cm(k)=dx(k)*opt5
+          dEvan1Calp(k)=0.0d0
+        enddo
+        Evan2=-opt6*Irsixp
+        do k=1,3
+          dEvan2Cat(k)=dx(k)*opt7
+          dEvan2Cm(k)=-dx(k)*opt7
+          dEvan2Calp(k)=0.0d0
+        enddo
+        ecation_prot=ecation_prot+Ec+Edip+Equad1+Equad2+Evan1+Evan2
+!        print *,ecation_prot,Ec+Edip+Equad1+Equad2+Evan1+Evan2
+        
+        do k=1,3
+          dEtotalCat(k)=dEcCat(k)+dEdipCat(k)+dEquad1Cat(k)+ &
+                       dEquad2Cat(k)+dEvan1Cat(k)+dEvan2Cat(k)
+!c             write(*,*) 'dEtotalCat inside', (dEtotalCat(l),l=1,3)
+          dEtotalCm(k)=dEcCm(k)+dEdipCm(k)+dEquad1Cm(k)+ &
+                      dEquad2Cm(k)+dEvan1Cm(k)+dEvan2Cm(k)
+          dEtotalCalp(k)=dEcCalp(k)+dEdipCalp(k)+dEquad1Calp(k) &
+                        +dEquad2Calp(k)+dEvan1Calp(k)+dEvan2Calp(k)
+        enddo
+            dscmag = 0.0d0
+            do k=1,3
+              dscvec(k) = dc(k,i+nres)
+              dscmag = dscmag+dscvec(k)*dscvec(k)
+            enddo
+            dscmag3 = dscmag
+            dscmag = sqrt(dscmag)
+            dscmag3 = dscmag3*dscmag
+            constA = 1.0d0+dASGL/dscmag
+            constB = 0.0d0
+            do k=1,3
+              constB = constB+dscvec(k)*dEtotalCm(k)
+            enddo
+            constB = constB*dASGL/dscmag3
+            do k=1,3
+              gg(k) = dEtotalCm(k)+dEtotalCalp(k)
+              gradpepcatx(k,i)=gradpepcatx(k,i)+ &
+               constA*dEtotalCm(k)-constB*dscvec(k)
+!            print *,j,constA,dEtotalCm(k),constB,dscvec(k)
+              gradpepcat(k,i)=gradpepcat(k,i)+gg(k)
+              gradpepcat(k,j)=gradpepcat(k,j)+dEtotalCat(k)
+             enddo
+        else if (itype(i,1).eq.13.or.itype(i,1).eq.14) then
+           if(itype(i,1).eq.14) then
+            inum=3
+            else
+            inum=4
+            endif
+            do k=1,6
+            vcatprm(k)=catprm(k,inum)
+            enddo
+            dASGL=catprm(7,inum)
+             do k=1,3
+                vcm(k)=(cm1(k)/cm1mag)*dASGL+c(k,i+nres)
+                valpha(k)=c(k,i)
+                vcat(k)=c(k,j)
+              enddo
 
-
+        do k=1,3
+          dx(k) = vcat(k)-vcm(k)
+        enddo
+        do k=1,3
+          v1(k)=(vcm(k)-valpha(k))
+          v2(k)=(vcat(k)-valpha(k))
+        enddo
+        v1m = sqrt(v1(1)**2+v1(2)**2+v1(3)**2)
+        v2m = sqrt(v2(1)**2+v2(2)**2+v2(3)**2)
+        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
+        wdip =vcatprm(2)
+        wdip=wdip/wh2o
+        wquad1 =vcatprm(3)
+        wquad1=wquad1/wh2o
+        wquad2 = vcatprm(4)
+        wquad2=wquad2/wh2o
+        wquad2p = 1-wquad2
+        wvan1 = vcatprm(5)
+        wvan2 =vcatprm(6)
+        opt = dx(1)**2+dx(2)**2
+        rsecp = opt+dx(3)**2
+        rs = sqrt(rsecp)
+        rthrp = rsecp*rs
+        rfourp = rthrp*rs
+        rsixp = rfourp*rsecp
+        reight=rsixp*rsecp
+        Ir = 1.0d0/rs
+        Irsecp = 1/rsecp
+        Irthrp = Irsecp/rs
+        Irfourp = Irthrp/rs
+        Irsixp = 1/rsixp
+        Ireight=1/reight
+        Irtw=Irsixp*Irsixp
+        Irthir=Irtw/rs
+        Irfourt=Irthir/rs
+        opt1 = (4*rs*dx(3)*wdip)
+        opt2 = 6*rsecp*wquad1*opt
+        opt3 = wquad1*wquad2p*Irsixp
+        opt4 = (wvan1*wvan2**12)
+        opt5 = opt4*12*Irfourt
+        opt6 = 2*wvan1*wvan2**6
+        opt7 = 6*opt6*Ireight
+        opt8 = wdip/v1m
+        opt10 = wdip/v2m
+        opt11 = (rsecp*v2m)**2
+        opt12 = (rsecp*v1m)**2
+        opt14 = (v1m*v2m*rsecp)**2
+        opt15 = -wquad1/v2m**2
+        opt16 = (rthrp*(v1m*v2m)**2)**2
+        opt17 = (v1m**2*rthrp)**2
+        opt18 = -wquad1/rthrp
+        opt19 = (v1m**2*v2m**2)**2
+        Edip=opt8*(v1dpv2)/(rsecp*v2m)
+        do k=1,3
+          dEdipCat(k)=opt8*(v1(k)*rsecp*v2m-((v2(k)/v2m&
+                     *rsecp+2*dx(k)*v2m)*v1dpv2))/opt11
+         dEdipCm(k)=opt10*(v2(k)*rsecp*v1m-((v1(k)/v1m&
+                    *rsecp-2*dx(k)*v1m)*v1dpv2))/opt12
+          dEdipCalp(k)=wdip*((-v1(k)-v2(k))*rsecp*v1m&
+                      *v2m-(-v1(k)/v1m*v2m*rsecp-v2(k)/v2m*v1m*rsecp)&
+                      *v1dpv2)/opt14
+        enddo
+        Equad1=-wquad1*v1dpv2**2/(rthrp*(v1m*v2m)**2)
+        do k=1,3
+          dEquad1Cat(k)=-wquad1*(2*v1(k)*v1dpv2*(rthrp*&
+                       (v1m*v2m)**2)-(3*dx(k)*rs*(v1m*v2m)**2+2*v1m*2*&
+                       v2(k)*1/2*1/v2m*v1m*v2m*rthrp)*v1dpv2**2)/opt16
+          dEquad1Cm(k)=-wquad1*(2*v2(k)*v1dpv2*(rthrp*&
+                      (v1m*v2m)**2)-(-3*dx(k)*rs*(v1m*v2m)**2+2*v2m*2*&
+                       v1(k)*1/2*1/v1m*v2m*v1m*rthrp)*v1dpv2**2)/opt16
+          dEquad1Calp(k)=opt18*(2*(-v1(k)-v2(k))*v1dpv2* &
+                        v1m**2*v2m**2-(-2*v1(k)*v2m**2-2*v2(k)*v1m**2)*&
+                        v1dpv2**2)/opt19
+        enddo
+        Equad2=wquad1*wquad2p*Irthrp
+        do k=1,3
+          dEquad2Cat(k)=-3*dx(k)*rs*opt3
+          dEquad2Cm(k)=3*dx(k)*rs*opt3
+          dEquad2Calp(k)=0.0d0
+        enddo
+        Evan1=opt4*Irtw
+        do k=1,3
+          dEvan1Cat(k)=-dx(k)*opt5
+          dEvan1Cm(k)=dx(k)*opt5
+          dEvan1Calp(k)=0.0d0
+        enddo
+        Evan2=-opt6*Irsixp
+        do k=1,3
+          dEvan2Cat(k)=dx(k)*opt7
+          dEvan2Cm(k)=-dx(k)*opt7
+          dEvan2Calp(k)=0.0d0
+        enddo
+         ecation_prot = ecation_prot+Edip+Equad1+Equad2+Evan1+Evan2
+        do k=1,3
+          dEtotalCat(k)=dEdipCat(k)+dEquad1Cat(k)+ &
+                       dEquad2Cat(k)+dEvan1Cat(k)+dEvan2Cat(k)
+          dEtotalCm(k)=dEdipCm(k)+dEquad1Cm(k)+ &
+                      dEquad2Cm(k)+dEvan1Cm(k)+dEvan2Cm(k)
+          dEtotalCalp(k)=dEdipCalp(k)+dEquad1Calp(k) &
+                        +dEquad2Calp(k)+dEvan1Calp(k)+dEvan2Calp(k)
+        enddo
+            dscmag = 0.0d0
+            do k=1,3
+              dscvec(k) = c(k,i+nres)-c(k,i)
+              dscmag = dscmag+dscvec(k)*dscvec(k)
+            enddo
+            dscmag3 = dscmag
+            dscmag = sqrt(dscmag)
+            dscmag3 = dscmag3*dscmag
+            constA = 1+dASGL/dscmag
+            constB = 0.0d0
+            do k=1,3
+              constB = constB+dscvec(k)*dEtotalCm(k)
+            enddo
+            constB = constB*dASGL/dscmag3
+            do k=1,3
+              gg(k) = dEtotalCm(k)+dEtotalCalp(k)
+              gradpepcatx(k,i)=gradpepcatx(k,i)+ &
+               constA*dEtotalCm(k)-constB*dscvec(k)
+              gradpepcat(k,i)=gradpepcat(k,i)+gg(k)
+              gradpepcat(k,j)=gradpepcat(k,j)+dEtotalCat(k)
+             enddo
+           else
+            rcal = 0.0d0
+            do k=1,3
+              r(k) = c(k,j)-c(k,i+nres)
+              rcal = rcal+r(k)*r(k)
+            enddo
+            ract=sqrt(rcal)
+            rocal=1.5
+            epscalc=0.2
+            r0p=0.5*(rocal+sig0(itype(i,1)))
+            r06 = r0p**6
+            r012 = r06*r06
+            Evan1=epscalc*(r012/rcal**6)
+            Evan2=epscalc*2*(r06/rcal**3)
+            r4 = rcal**4
+            r7 = rcal**7
+            do k=1,3
+              dEvan1Cm(k) = 12*r(k)*epscalc*r012/r7
+              dEvan2Cm(k) = 12*r(k)*epscalc*r06/r4
+            enddo
+            do k=1,3
+              dEtotalCm(k)=dEvan1Cm(k)+dEvan2Cm(k)
+            enddo
+                 ecation_prot = ecation_prot+ Evan1+Evan2
+            do  k=1,3
+               gradpepcatx(k,i)=gradpepcatx(k,i)+ & 
+               dEtotalCm(k)
+              gradpepcat(k,i)=gradpepcat(k,i)+dEtotalCm(k)
+              gradpepcat(k,j)=gradpepcat(k,j)-dEtotalCm(k)
+             enddo
+         endif ! 13-16 residues
+       enddo !j
+       enddo !i
        return
        end subroutine ecat_prot
 

diff --git a/source/unres/io.f90 b/source/unres/io.f90
index a86fc0b..dbc5ce0 100644 (file)
--- a/source/unres/io.f90
+++ b/source/unres/io.f90
@@ -798,6 +798,8 @@
       call reada(weightcard,'TEMP0',temp0,300.0d0)
       if (index(weightcard,'SOFT').gt.0) ipot=6
       call reada(weightcard,'WBOND_NUCL',wbond_nucl,1.0D0)
+      call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
+      call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
 ! 12/1/95 Added weight for the multi-body term WCORR
       call reada(weightcard,'WCORRH',wcorr,1.0D0)
       if (wcorr4.gt.0.0d0) wcorr=wcorr4

diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90
index f6a3919..8e913d0 100644 (file)
--- a/source/unres/io_config.f90
+++ b/source/unres/io_config.f90
@@ -739,7 +739,7 @@
       character(len=3) :: lancuch      !,ucase
 !el  local variables
       integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
-      integer :: maxinter,junk,kk,ii
+      integer :: maxinter,junk,kk,ii,ncatprotparm
       real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
                 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
                 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
@@ -854,6 +854,19 @@
           enddo
         enddo
       endif
+            do i=1,ntyp_molec(5)
+             read(iion,*) msc(i,5),restok(i,5)
+             print *,msc(i,5),restok(i,5)
+            enddo
+            ip(5)=0.2
+!            isc(5)=0.2
+            read (iion,*) ncatprotparm
+            allocate(catprm(ncatprotparm,4))
+            do k=1,4
+            read (iion,*)  (catprm(i,k),i=1,ncatprotparm)
+            enddo
+            print *, catprm
+!            read (iion,*) (vcatprm(k),k=1,ncatprotpram)
 !----------------------------------------------------
       allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
       allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp))     !(-ntyp:ntyp)
@@ -2578,7 +2591,7 @@
       character(len=80) :: card
       real(kind=8),dimension(3,20) :: sccor
       integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin      !rescode,
-      integer :: isugar,molecprev
+      integer :: isugar,molecprev,firstion
       character*1 :: sugar
       real(kind=8) :: cou
       real(kind=8),dimension(3) :: ccc
@@ -2715,6 +2728,7 @@
              if (atom.eq.'CA  '.or.atom.eq.'N   ') then
              molecule=1
               itype(ires,molecule)=rescode(ires,res,0,molecule)
+              firstion=0
 !              nres_molec(molecule)=nres_molec(molecule)+1
             else
              molecule=2
@@ -2771,6 +2785,7 @@
              if (counter.eq.5) then
 !            iii=iii+1
               nres_molec(molecule)=nres_molec(molecule)+1
+              firstion=0
 !            do j=1,3
 !              sccor(j,iii)=c(j,ires)
 !            enddo
@@ -2810,7 +2825,16 @@
               endif
           endif
         else if ((ions).and.(card(1:6).eq.'HETATM')) then
-         
+       if (firstion.eq.0) then 
+       firstion=1
+       if (unres_pdb) then
+         do j=1,3
+           dc(j,ires)=sccor(j,iii)
+         enddo
+       else
+          call sccenter(ires,iii,sccor)
+       endif
+       endif
           read (card(12:16),*) atom
           print *,"HETATOM", atom
           read (card(18:20),'(a3)') res
@@ -2822,6 +2846,7 @@
            if (molecule.ne.5) molecprev=molecule
            molecule=5
            nres_molec(molecule)=nres_molec(molecule)+1
+           print *,"HERE",nres_molec(molecule)
            itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
            read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
           endif
@@ -2832,7 +2857,8 @@
 ! Calculate dummy residue coordinates inside the "chain" of a multichain
 ! system
       nres=ires
-      if (ires_old.ne.ires) nres_molec(molecule)=nres_molec(molecule)-2
+      if ((ires_old.ne.ires).and.(molecule.ne.5)) &
+         nres_molec(molecule)=nres_molec(molecule)-2
 !      print *,'I have', nres_molec(:)
       
       do k=1,4 ! ions are without dummy 
@@ -3084,6 +3110,7 @@
 
           enddo
           itype_temporary(seqalingbegin,k)=itype(i,k)
+          print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
           istype_temp(seqalingbegin)=istype(i)
           molnum(seqalingbegin)=k
           seqalingbegin=seqalingbegin+1
@@ -4284,6 +4311,8 @@
       open (iliptranpar,file=liptranname,status='old',action='read')
       call getenv_loc('TUBEPAR',tubename)
       open (itube,file=tubename,status='old',action='read')
+      call getenv_loc('IONPAR',ionname)
+      open (iion,file=ionname,status='old',action='read')
 
 !      print *,"Processor",myrank," opened file ISIDEP" 
 !      print *,"Processor",myrank," opened parameter files" 
@@ -4328,6 +4357,8 @@
       open (iliptranpar,file=liptranname,status='old')
       call getenv_loc('TUBEPAR',tubename)
       open (itube,file=tubename,status='old')
+      call getenv_loc('IONPAR',ionname)
+      open (iion,file=ionname,status='old')
 #else
       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
         readonly)
@@ -4376,6 +4407,8 @@
       open (iliptranpar,file=liptranname,status='old',action='read')
       call getenv_loc('TUBEPAR',tubename)
       open (itube,file=tubename,status='old',action='read')
+      call getenv_loc('IONPAR',ionname)
+      open (iion,file=ionname,status='old',action='read')
 
 #ifndef CRYST_SC
       call getenv_loc('ROTPARPDB',rotname_pdb)