working cutoff
authorAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Wed, 18 Mar 2020 11:52:31 +0000 (12:52 +0100)
committerAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Wed, 18 Mar 2020 11:52:31 +0000 (12:52 +0100)
source/unres/CMakeLists.txt
source/unres/control.F90
source/unres/data/energy_data.F90
source/unres/energy.F90
source/unres/io_config.F90
source/wham/conform_compar.F90
source/wham/io_database.F90
source/wham/io_wham.F90
source/wham/wham_data.F90

index 0d12454..3b6f103 100644 (file)
@@ -72,7 +72,7 @@ if (Fortran_COMPILER_NAME STREQUAL "ifort")
   set (CMAKE_Fortran_FLAGS_RELEASE " ")
   set (CMAKE_Fortran_FLAGS_DEBUG   "-O0 -g -traceback")
 #  set(FFLAGS0 "-fpp -c -O3 -ip " ) 
-  set(FFLAGS0 "-CB -g -ip -fpp  -heap-arrays" ) 
+  set(FFLAGS0 "-O3 -ip -fpp  -heap-arrays -mcmodel=medium" ) 
 #  set(FFLAGS0 "-O0 -CB -CA -g" )
   set(FFLAGS1 "-fpp -c -O " ) 
   set(FFLAGS2 "-fpp -c -g -CA -CB ")
index a1e49d3..b46dd40 100644 (file)
 !
 ! The following is just to define auxiliary variables used in angle conversion
 !
+!      ifirstrun=0
       pi=4.0D0*datan(1.0D0)
       dwapi=2.0D0*pi
       dwapi3=dwapi/3.0D0
index b0995af..8fd85cd 100644 (file)
         real(kind=8),dimension(:,:,:),allocatable :: alphasurcat,&
            alphisocat,wqdipcat,dtailcat,wstatecat
          real(kind=8),dimension(:,:,:,:),allocatable :: dheadcat
-
+          integer,dimension(60000) :: contlistscpi_f,contlistscpj_f
+!         integer :: ifirstrun,ilist_scp_first
 !        real(kind=8),dimension(:,:),allocatable :: alphapol,epshead,&
 !           sig0head,sigiso1,sigiso2,rborn,sigmap1,sigmap2,chis,wquad,chipp,&
 !           epsintab,debaykap
 
 
 !end of ions parameters by Agnieszka Lipska (Ca, K, Na, Mg, Cl)-----------------------
-
+!
+! FRAGMENT FOR INTERACTION LIST
+        integer,dimension(:),allocatable :: newcontlistppi,newcontlistppj,&
+        newcontlisti,newcontlistj,  newcontlistscpi,newcontlistscpj
+        integer :: g_listpp_start,g_listpp_end,g_listscp_start,g_listscp_end,&
+        g_listscsc_start,g_listscsc_end   
       end module energy_data
index 8d39a59..b0a94ad 100644 (file)
        integer ishield_listbuf(-1:nres), &
        shield_listbuf(maxcontsshi,-1:nres),k,j,i,iii,impishi,mojint,jjj
 !       print *,"I START ENERGY"
-       imatupdate=1
+       imatupdate=100
 !       if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
 !      real(kind=8),  dimension(:),allocatable::  fac_shieldbuf 
 !      real(kind=8), dimension(:,:,:),allocatable:: &
         time_Bcastw=time_Bcastw+MPI_Wtime()-time00
 !        call chainbuild_cart
       endif
-!       if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
+!       print *,"itime_mat",itime_mat,imatupdate
+        if (nfgtasks.gt.1) then 
+        call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR)
+        endif
+       if (mod(itime_mat,imatupdate).eq.0) call make_SCp_inter_list
+       if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
+       if (mod(itime_mat,imatupdate).eq.0) call make_pp_inter_list
 
 !      print *,'Processor',myrank,' calling etotal ipot=',ipot
 !      print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
        etors_nucl=0.0d0
        estr_nucl=0.0d0
        ecorr3_nucl=0.0d0
+       ecorr_nucl=0.0d0
        ebe_nucl=0.0d0
        evdwsb=0.0d0
        eelsb=0.0d0
        eelpsb=0.0d0
        evdwpp=0.0d0
        eespp=0.0d0
+       etors_d_nucl=0.0d0
       endif
 !      write(iout,*) ecorr_nucl,"ecorr_nucl",nres_molec(2)
 !      print *,"before ecatcat",wcatcat
+      if (nres_molec(5).gt.0) then
       if (nfgtasks.gt.1) then
       if (fg_rank.eq.0) then
       call ecatcat(ecationcation)
       else
       call ecats_prot_amber(ecation_prot)
       endif
-      if (nres_molec(2).gt.0) then
+      else
+      ecationcation=0.0d0
+      ecation_prot=0.0d0
+      endif
+      if ((nres_molec(2).gt.0).and.(nres_molec(1).gt.0)) then
       call eprot_sc_base(escbase)
       call epep_sc_base(epepbase)
       call eprot_sc_phosphate(escpho)
 !      include 'COMMON.SBRIDGE'
       logical :: lprn
 !el local variables
-      integer :: iint,itypi,itypi1,itypj,subchap
+      integer :: iint,itypi,itypi1,itypj,subchap,icont
       real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
       real(kind=8) :: evdw,sig0ij
       real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
       dPOLdOM1=0.0d0
 
 
-      do i=iatsc_s,iatsc_e
+      do icont=g_listscsc_start,g_listscsc_end
+      i=newcontlisti(icont)
+      j=newcontlistj(icont)
+
+!      do i=iatsc_s,iatsc_e
 !C        print *,"I am in EVDW",i
         itypi=iabs(itype(i,1))
 !        if (i.ne.47) cycle
 !
 ! Calculate SC interaction energy.
 !
-        do iint=1,nint_gr(i)
-          do j=istart(i,iint),iend(i,iint)
+!        do iint=1,nint_gr(i)
+!          do j=istart(i,iint),iend(i,iint)
             IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
               call dyn_ssbond_ene(i,j,evdwij)
               evdw=evdw+evdwij
 ! Calculate angular part of the gradient.
             call sc_grad
             ENDIF    ! dyn_ss            
-          enddo      ! j
-        enddo        ! iint
+!          enddo      ! j
+!        enddo        ! iint
       enddo          ! i
 !       print *,"ZALAMKA", evdw
 !      write (iout,*) "Number of loop steps in EGB:",ind
                                              0.0d0,1.0d0,0.0d0,&
                                              0.0d0,0.0d0,1.0d0/),shape(unmat)) 
 !el local variables
-      integer :: i,k,j
+      integer :: i,k,j,icont
       real(kind=8) :: ees,evdw1,eel_loc,eello_turn3,eello_turn4
       real(kind=8) :: fac,t_eelecij,fracinbuf
     
 ! Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 !
 !      print *,"iatel_s,iatel_e,",iatel_s,iatel_e
-      do i=iatel_s,iatel_e
+!      do i=iatel_s,iatel_e
+! JPRDLC
+       do icont=g_listpp_start,g_listpp_end
+        i=newcontlistppi(icont)
+        j=newcontlistppj(icont)
         if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
 
 !        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         num_conti=num_cont_hb(i)
-        do j=ielstart(i),ielend(i)
+!        do j=ielstart(i),ielend(i)
 !          write (iout,*) i,j,itype(i,1),itype(j,1)
           if (itype(j,1).eq.ntyp1.or. itype(j+1,1).eq.ntyp1) cycle
           call eelecij(i,j,ees,evdw1,eel_loc)
-        enddo ! j
+!        enddo ! j
         num_cont_hb(i)=num_conti
       enddo   ! i
 !      write (iout,*) "Number of loop steps in EELEC:",ind
 !      include 'COMMON.CONTROL'
       real(kind=8),dimension(3) :: ggg
 !el local variables
-      integer :: i,iint,j,k,iteli,itypj,subchap
+      integer :: i,iint,j,k,iteli,itypj,subchap,icont
       real(kind=8) :: evdw2,evdw2_14,xi,yi,zi,xj,yj,zj,rrij,fac,&
                    e1,e2,evdwij,rij
       real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
       evdw2_14=0.0d0
 !d    print '(a)','Enter ESCP'
 !d    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
-      do i=iatscp_s,iatscp_e
+!      do i=iatscp_s,iatscp_e
+       do icont=g_listscp_start,g_listscp_end
+        i=newcontlistscpi(icont)
+        j=newcontlistscpj(icont)
         if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
         iteli=itel(i)
         xi=0.5D0*(c(1,i)+c(1,i+1))
           zi=mod(zi,boxzsize)
           if (zi.lt.0) zi=zi+boxzsize
 
-        do iint=1,nscp_gr(i)
+!        do iint=1,nscp_gr(i)
 
-        do j=iscpstart(i,iint),iscpend(i,iint)
+!        do j=iscpstart(i,iint),iscpend(i,iint)
           itypj=iabs(itype(j,1))
           if (itypj.eq.ntyp1) cycle
 ! Uncomment following three lines for SC-p interactions
             gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
             gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
           enddo
-        enddo
+!        enddo
 
-        enddo ! iint
+!        enddo ! iint
       enddo ! i
       do i=1,nct
         do j=1,3
       allocate(uygrad(3,3,2,nres))
       allocate(uzgrad(3,3,2,nres))
 !(3,3,2,maxres)
+! allocateion of lists JPRDLA
+      allocate(newcontlistppi(200*nres))
+      allocate(newcontlistscpi(200*nres))
+      allocate(newcontlisti(200*nres))
+      allocate(newcontlistppj(200*nres))
+      allocate(newcontlistscpj(200*nres))
+      allocate(newcontlistj(200*nres))
 
       return
       end subroutine alloc_ener_arrays
       include 'mpif.h'
       real*8 :: xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp
       real*8 :: dist_init, dist_temp,r_buff_list
-      integer:: contlisti(50*nres),contlistj(50*nres)
-      integer :: newcontlisti(50*nres),newcontlistj(50*nres) 
+      integer:: contlisti(200*nres),contlistj(200*nres)
+!      integer :: newcontlisti(200*nres),newcontlistj(200*nres) 
       integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_sc,g_ilist_sc
       integer displ(0:nprocs),i_ilist_sc(0:nprocs),ierr
 !            print *,"START make_SC"
+          r_buff_list=5.0
             ilist_sc=0
             do i=iatsc_s,iatsc_e
              itypi=iabs(itype(i,1))
       write (iout,*) i,newcontlisti(i),newcontlistj(i)
       enddo
 #endif
-
+        call int_bounds(g_ilist_sc,g_listscsc_start,g_listscsc_end)
       return
       end subroutine make_SCSC_inter_list
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
       subroutine make_SCp_inter_list
+      use MD_data,  only: itime_mat
+
       include 'mpif.h'
       real*8 :: xi,yi,zi,xj,yj,zj,xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp
       real*8 :: dist_init, dist_temp,r_buff_list
-      integer:: contlistscpi(50*nres),contlistscpj(50*nres)
-      integer :: newcontlistscpi(50*nres),newcontlistscpj(50*nres)
+      integer:: contlistscpi(200*nres),contlistscpj(200*nres)
+!      integer :: newcontlistscpi(200*nres),newcontlistscpj(200*nres)
       integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_scp,g_ilist_scp
       integer displ(0:nprocs),i_ilist_scp(0:nprocs),ierr
-            print *,"START make_SC"
+!            print *,"START make_SC"
+      r_buff_list=5.0
             ilist_scp=0
       do i=iatscp_s,iatscp_e
         if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
           yj=yj_safe-yi
           zj=zj_safe-zi
        endif
+#ifdef DEBUG
+                ! r_buff_list is a read value for a buffer 
+               if ((sqrt(dist_init).le.(r_cut_ele)).and.(ifirstrun.eq.0)) then
+! Here the list is created
+                 ilist_scp_first=ilist_scp_first+1
+! this can be substituted by cantor and anti-cantor
+                 contlistscpi_f(ilist_scp_first)=i
+                 contlistscpj_f(ilist_scp_first)=j
+              endif
+#endif
 ! r_buff_list is a read value for a buffer 
                if (sqrt(dist_init).le.(r_cut_ele+r_buff_list)) then
 ! Here the list is created
       do i=1,g_ilist_scp
       write (iout,*) i,newcontlistscpi(i),newcontlistscpj(i)
       enddo
+
+!      if (ifirstrun.eq.0) ifirstrun=1
+!      do i=1,ilist_scp_first
+!       do j=1,g_ilist_scp
+!        if ((newcontlistscpi(j).eq.contlistscpi_f(i)).and.&
+!         (newcontlistscpj(j).eq.contlistscpj_f(i))) go to 126
+!        enddo
+!       print *,itime_mat,"ERROR matrix needs updating"
+!       print *,contlistscpi_f(i),contlistscpj_f(i)
+!  126  continue
+!      enddo
 #endif
+        call int_bounds(g_ilist_scp,g_listscp_start,g_listscp_end)
+
       return
       end subroutine make_SCp_inter_list
 
       real*8 :: xmedj,ymedj,zmedj
       real*8 :: dist_init, dist_temp,r_buff_list,dxi,dyi,dzi,xmedi,ymedi,zmedi
       real*8 :: dx_normi,dy_normi,dz_normi,dxj,dyj,dzj,dx_normj,dy_normj,dz_normj
-      integer:: contlistppi(50*nres),contlistppj(50*nres)
-      integer :: newcontlistppi(50*nres),newcontlistppj(50*nres)
+      integer:: contlistppi(200*nres),contlistppj(200*nres)
+!      integer :: newcontlistppi(200*nres),newcontlistppj(200*nres)
       integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_pp,g_ilist_pp
       integer displ(0:nprocs),i_ilist_pp(0:nprocs),ierr
-            print *,"START make_SC"
+!            print *,"START make_SC"
             ilist_pp=0
+      r_buff_list=5.0
       do i=iatel_s,iatel_e
         if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
         dxi=dc(1,i)
         newcontlistppj(i)=contlistppj(i)
         enddo
         endif
-
+        call int_bounds(g_ilist_pp,g_listpp_start,g_listpp_end)
 #ifdef DEBUG
       write (iout,*) "after MPIREDUCE",g_ilist_pp
       do i=1,g_ilist_pp
index 847df40..c37e8a5 100644 (file)
       endif
 
 ! CUTOFFF ON ELECTROSTATICS
-      call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
+      call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
       write(iout,*) "R_CUT_ELE=",r_cut_ele
 ! Lipidic parameters
index 1720af2..1d3a632 100644 (file)
@@ -25,7 +25,7 @@
 !-----------------------------------------------------------------------------
 ! conf_compar.F
 !-----------------------------------------------------------------------------
-      subroutine conf_compar(jcon,lprn,print_class)
+      subroutine conf_compar(jcon)
 !      implicit real*8 (a-h,o-z)
       use energy_data, only:icont,ncont,nnt,nct,maxcont!,&
 !                      nsccont_frag_ref,isccont_frag_ref
                  nc_match,ncon_match,iclass_rms,ishifft_rms,ishiff,ishif
       integer :: k,kk,iclass_con,iscor,ik,ishifft_con,idig,iex,im
 !      print *,"Enter conf_compar",jcon
-      call angnorm12(rmsang)
+       lprn=.false.
+       print_class=.false.
+        write (iout,*) "before anything"
+        call flush(iout)
+
+!      call angnorm12(rmsang)
 ! Level 1: check secondary and supersecondary structure
-      call elecont(lprn,ncont,icont,nnt,nct)
+!      call elecont(lprn,ncont,icont,nnt,nct)
       if (lprn) then
         write (iout,*) "elecont finished"
         call flush(iout)
         ymedi=yi+0.5*dyi
         zmedi=zi+0.5*dzi
         do 4 j=i+2,ien-1
-          ind=ind+1
+!          ind=ind+1
           iteli=itel(i)
           itelj=itel(j)
           if (j.eq.i+2 .and. itelj.eq.2) iteli=2
index 81cd3f3..5b8cad4 100644 (file)
@@ -442,7 +442,7 @@ write(iout,*) "end of read database"
 !      rtime=0.0d0
 !      rpotE=0.0d0
 !      rt_bath=0.0d0
-
+      rmsdev=0.0d0
       call set_slices(is,ie,ts,te,iR,ib,iparm)
       nprop_prev=0
       do i=1,nQ
@@ -916,7 +916,8 @@ write(iout,*) "end of read database"
       use geometry_data
       use control_data, only:indpdb
       use w_compar_data
-      use conform_compar, only:conf_compar
+      use conform_compar, only:conf_compar,rmsnat,qwolynes
+      use energy_data, only:icont,ncont,nnt,nct,maxcont!,&
 !      implicit none
 !      include "DIMENSIONS"
 !      include "DIMENSIONS.ZSCOPT"
@@ -954,13 +955,15 @@ write(iout,*) "end of read database"
       integer :: i,itj,ii,iii,j,k,l
       integer :: ixdrf,iret
       integer :: iscor,islice
-      real(kind=8) :: rmsdev,efree,eini
+      real(kind=8) :: rmsdev,efree,eini,qnat2
       real(kind=4) :: csingle(3,nres*2)
       real(kind=8) :: energ
+       
 !      integer ilen,iroof
 !      external ilen,iroof
       integer :: ir,ib,iparm
       integer :: isecstr(nres)
+      logical :: test
       write (licz2,'(bz,i2.2)') islice
       call opentmp(islice,ientout,bprotfile_temp)
       write (iout,*) "bprotfile_temp ",bprotfile_temp
@@ -1041,12 +1044,14 @@ write(iout,*) "end of read database"
         print *,istat,statname
         open(istat,file=statname,status="unknown")
       endif
-
+      print *,"Tu dochodze"
+      print *,scount(me)
 #ifdef MPI
       do i=1,scount(me)
 #else
       do i=1,ntot(islice)
 #endif
+        print *,"before ientout read"
         read(ientout,rec=i,err=101) &
           ((csingle(l,k),l=1,3),k=1,nres),&
           ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
@@ -1063,14 +1068,24 @@ write(iout,*) "end of read database"
 !        write (iout,*) "Calling conf_compar",i
 !        call flush(iout)
          anatemp= 1.0d0/(beta_h(ib,iparm)*1.987D-3)
+        print *,"before conf_compar"
         if (indpdb.gt.0) then
-          call conf_compar(i,.false.,.true.)
+        print *,"just before conf_compar",i
+        print *,icont,ncont,nnt,nct,"maxcont",maxcont
+        test=.false.
+!          call conf_compar(i,.false.,.true.)
+!          call conf_compar(i)
+!           call rmsnat(i)
+           rms_nat=rmsnat(i)
+           qnat2=qwolynes(0,0) 
+         print *,"just after conf_compar"
 !        else
 !            call elecont(.false.,ncont,icont,nnt,nct)
 !            call secondary2(.false.,.false.,ncont,icont,isecstr)
         endif
 !        write (iout,*) "Exit conf_compar",i
 !        call flush(iout)
+         print *,"before ientin"
         if (bxfile .or.cxfile .or. ensembles.gt.0) write (ientin,rec=i) &
           ((csingle(l,k),l=1,3),k=1,nres),&
           ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
@@ -1087,6 +1102,7 @@ write(iout,*) "end of read database"
       close(istat)
       if (bxfile .or. cxfile .or. ensembles.gt.0) close(ientin)
 #ifdef MPI
+      print *,"before MPI_barrier"
       call MPI_Barrier(WHAM_COMM,IERROR)
       if (me.ne.Master .or. .not.bxfile .and. .not. cxfile &
          .and. ensembles.eq.0) return
@@ -1307,7 +1323,7 @@ write(iout,*) "end of read database"
       integer :: islice,iR,ib,iparm
       integer :: is(MaxSlice),ie(MaxSlice),nrec_slice
       real(kind=8) :: ts(MaxSlice),te(MaxSlice),time_slice
-
+      time_slice=0
       do islice=1,nslice
         if (time_end_collect(iR,ib,iparm).ge.1.0d10) then
           ts(islice)=time_start_collect(iR,ib,iparm)
index cb2e1ec..b96b816 100644 (file)
@@ -667,6 +667,7 @@ allocate(ww(max_eneW))
       endif
             if (.not. allocated(msc)) allocate(msc(ntyp1,5))
             if (.not. allocated(restok)) allocate(restok(ntyp1,5))
+       if (oldion.eq.1) then
 
             do i=1,ntyp_molec(5)
              read(iion,*) msc(i,5),restok(i,5)
@@ -680,6 +681,7 @@ allocate(ww(max_eneW))
             read (iion,*)  (catprm(i,k),i=1,ncatprotparm)
             enddo
             print *, catprm
+      endif
       read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
       do i=1,ntyp_molec(2)
         nbondterm_nucl(i)=1
@@ -2671,6 +2673,125 @@ allocate(ww(max_eneW))
           bad(i,j)=0.0D0
         enddo
       enddo
+! Ions by Aga
+
+       allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
+       allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
+       allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
+       allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
+       allocate(epsintabcat(ntyp,ntyp))
+       allocate(dtailcat(2,ntyp,ntyp))
+       allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
+       allocate(wqdipcat(2,ntyp,ntyp))
+       allocate(wstatecat(4,ntyp,ntyp))
+       allocate(dheadcat(2,2,ntyp,ntyp))
+       allocate(nstatecat(ntyp,ntyp))
+       allocate(debaykapcat(ntyp,ntyp))
+
+      if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
+      if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
+!      if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
+      if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
+      if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
+
+
+      allocate (ichargecat(ntyp_molec(5)))
+! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
+       if (oldion.eq.0) then
+            if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
+            allocate(icharge(1:ntyp1))
+            read(iion,*) (icharge(i),i=1,ntyp)
+            else
+             read(iion,*) ijunk
+            endif
+
+            do i=1,ntyp_molec(5)
+             read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
+             print *,msc(i,5),restok(i,5)
+            enddo
+            ip(5)=0.2
+!DIR$ NOUNROLL 
+      do j=1,ntyp_molec(5)
+       do i=1,ntyp
+!       do j=1,ntyp_molec(5)
+!        write (*,*) "Im in ALAB", i, " ", j
+        read(iion,*) &
+       epscat(i,j),sigmacat(i,j), &
+!       chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
+       chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
+
+       (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
+!       chiscat(i,j),chiscat(j,i), &
+       chis1cat(i,j),chis2cat(i,j), &
+
+       nstatecat(i,j),(wstatecat(k,i,j),k=1,4), &                           !5 w tej lini - 1 integer pierwszy
+       dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
+       dtailcat(1,i,j),dtailcat(2,i,j), &
+       epsheadcat(i,j),sig0headcat(i,j), &
+!wdipcat = w1 , w2
+!       rborncat(i,j),rborncat(j,i),&
+       rborn1cat(i,j),rborn2cat(i,j),&
+       (wqdipcat(k,i,j),k=1,2), &
+       alphapolcat(i,j),alphapolcat(j,i), &
+       (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
+!       print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i) 
+       END DO
+      END DO
+      allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
+      do i=1,ntyp
+        do j=1,ntyp_molec(5)
+          epsij=epscat(i,j)
+          rrij=sigmacat(i,j)
+          rrij=rrij**expon
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_aq_cat(i,j)=epsij*rrij*rrij
+          bb_aq_cat(i,j)=-sigeps*epsij*rrij
+         enddo
+       enddo
+       do i=1,ntyp
+       do j=1,ntyp_molec(5)
+      if (i.eq.10) then
+      write (iout,*) 'i= ', i, ' j= ', j
+      write (iout,*) 'epsi0= ', epscat(i,j)
+      write (iout,*) 'sigma0= ', sigmacat(i,j)
+      write (iout,*) 'chi1= ', chi1cat(i,j)
+      write (iout,*) 'chi1= ', chi2cat(i,j)
+      write (iout,*) 'chip1= ', chipp1cat(1,j)
+      write (iout,*) 'chip2= ', chipp2cat(1,j)
+      write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
+      write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
+      write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
+      write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
+      write (iout,*) 'sig1= ', sigmap1cat(1,j)
+      write (iout,*) 'sig2= ', sigmap2cat(1,j)
+      write (iout,*) 'chis1= ', chis1cat(1,j)
+      write (iout,*) 'chis1= ', chis2cat(1,j)
+      write (iout,*) 'nstatecat(i,j)= ', nstatecat(1,j)
+      write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,1,j)
+      write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
+      write (iout,*) 'dhead2= ', dheadcat(1,2,1,j)
+      write (iout,*) 'a1= ', rborn1cat(i,j)
+      write (iout,*) 'a2= ', rborn2cat(i,j)
+      write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
+      write (iout,*) 'alphapol1= ',  alphapolcat(1,j)
+      write (iout,*) 'alphapol2= ',  alphapolcat(j,1)
+      write (iout,*) 'w1= ', wqdipcat(1,i,j)
+      write (iout,*) 'w2= ', wqdipcat(2,i,j)
+      write (iout,*) 'debaykapcat(i,j)= ',  debaykapcat(1,j)
+      endif
+
+      If ((i.eq.1).and.(j.eq.27)) then
+      write (iout,*) 'SEP'
+      Write (iout,*) 'w1= ', wqdipcat(1,1,27)
+      Write (iout,*) 'w2= ', wqdipcat(2,1,27)
+      endif
+
+       enddo
+       enddo
+
+      endif
+
 #ifdef CLUSTER
 !
 ! Define the SC-p interaction constants
@@ -2974,6 +3095,8 @@ allocate(ww(max_eneW))
       call reada(controlcard,'BOXY',boxysize,100.0d0)
       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
       call readi(controlcard,"SCELEMODE",scelemode,0)
+      call readi(controlcard,"OLDION",oldion,0)
+
       print *,"SCELE",scelemode
       call readi(controlcard,'TORMODE',tor_mode,0)
 !C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
index f505680..28fcbda 100644 (file)
       logical :: punch_dist,print_rms,caonly,verbose,merge_helices,&
         bxfile,cxfile,histfile,entfile,zscfile,rmsrgymap,&
         with_dihed_constr,check_conf,histout
-      integer :: icomparfunc,pdbint,ensembles,constr_dist
+      integer :: icomparfunc,pdbint,ensembles,constr_dist,oldion
 !---------------------------------------------------------------------------
 ! COMMON.OBCINKA
 !      common /obcinka/