correction in new ions and clusterfix
[unres4.git] / source / unres / io_config.F90
index 345ec71..61cb7b4 100644 (file)
           enddo
         enddo
       endif
+
+
+
+      if (.not.allocated(ichargecat)) allocate (ichargecat(ntyp_molec(5)))
        if (oldion.eq.1) then
             do i=1,ntyp_molec(5)
-             read(iion,*) msc(i,5),restok(i,5)
+             read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
              print *,msc(i,5),restok(i,5)
             enddo
             ip(5)=0.2
             enddo
             print *, catprm
          endif
+      allocate(catnuclprm(14,ntyp_molec(2),ntyp_molec(5)))
+      do i=1,ntyp_molec(5)
+         do j=1,ntyp_molec(2)
+         write(iout,*) i,j
+            read(iionnucl,*) (catnuclprm(k,j,i),k=1,14)
+         enddo
+      enddo
+      write(*,'(3(5x,a6)11(7x,a6))') "w1    ","w2    ","epslj ","pis1  ", &
+      "sigma0","epsi0 ","chi1   ","chip1 ","sig   ","b1    ","b2    ", &
+      "b3    ","b4    ","chis1  "
+      do i=1,ntyp_molec(5)
+         do j=1,ntyp_molec(2)
+            write(*,'(3(f10.3,x),11(f12.6,x),a3,2a)') (catnuclprm(k,j,i),k=1,14), &
+                                      restyp(i,5),"-",restyp(j,2)
+         enddo
+      enddo
 !            read (iion,*) (vcatprm(k),k=1,ncatprotpram)
 !----------------------------------------------------
       allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
        (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &           !6 w tej linii
        chis(i,j),chis(j,i)                                            !2 w tej linii
         endif
-!       print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i) 
+       print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i), wqdip(1,i,j)
        END DO
       END DO
       do i=1,ntyp
 ! Ions by Aga
 
        allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
-       allocate(sigiso1cat(ntyp,ntyp),rborncat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
+       allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
        allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
-       allocate(chiscat(ntyp,ntyp),wquadcat(ntyp,ntyp),chippcat(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(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(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)
+
 
+            if (.not.allocated(ichargecat)) 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)
+             read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
              print *,msc(i,5),restok(i,5)
             enddo
             ip(5)=0.2
-
-      do i=1,ntyp
-       do j=1,ntyp_molec(5)
+!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), &
+       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), &
+!       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),(wqdipcat(k,i,j),k=1,2), &
+!       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) 
+!     if (i.eq.1) then
+!     write (iout,*) 'i= ', i, ' j= ', j
+!     write (iout,*) 'epsi0= ', epscat(1,j)
+!     write (iout,*) 'sigma0= ', sigmacat(1,j)
+!     write (iout,*) 'chi(i,j)= ', chicat(1,j)
+!     write (iout,*) 'chi(j,i)= ', chicat(j,1)
+!     write (iout,*) 'chip(i,j)= ', chippcat(1,j)
+!     write (iout,*) 'chip(j,i)= ', chippcat(j,1)
+!     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,*) 'chis(i,j)= ', chiscat(1,j)
+!     write (iout,*) 'chis(j,i)= ', chiscat(j,1)
+!     write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
+!     write (iout,*) 'a1= ', rborncat(j,1)
+!     write (iout,*) 'a2= ', rborncat(1,j)
+!     write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
+!     write (iout,*) 'alphapol1= ',  alphapolcat(1,j)
+!     write (iout,*) 'w1= ', wqdipcat(1,1,j)
+!     write (iout,*) 'w2= ', wqdipcat(2,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
+
        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
 
       
          istype(i)=istype_temp(i)
         enddo
        enddo
+       if ((itype(1,1).eq.ntyp1).and.itype(2,5).ne.0) then
+! I have only ions now dummy atoms in the system        
+       molnum(1)=5
+       itype(1,5)=itype(2,5)
+       itype(1,1)=0
+       do i=2,nres
+         itype(i,5)=itype(i+1,5)
+       enddo
+       itype(nres,5)=0
+       nres=nres-1
+       nres_molec(1)=nres_molec(1)-1
+      endif
 !      if (itype(1,1).eq.ntyp1) then
 !        nsup=nsup-1
 !        nstart_sup=2
       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
       protein=index(controlcard,"PROTEIN").gt.0
       ions=index(controlcard,"IONS").gt.0
+      call readi(controlcard,'OLDION',oldion,1)
       nucleic=index(controlcard,"NUCLEIC").gt.0
       write (iout,*) "with_theta_constr ",with_theta_constr
       AFMlog=(index(controlcard,'AFM'))
       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
 !      enddo
       buff_shield=1.0d0
       endif
+      itime_mat=0
       return
       end subroutine read_control
 !-----------------------------------------------------------------------------
       open (itube,file=tubename,status='old')
       call getenv_loc('IONPAR',ionname)
       open (iion,file=ionname,status='old')
+      call getenv_loc('IONPAR_NUCL',ionnuclname)
+      open (iionnucl,file=ionnuclname,status='old')
 #else
       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
         readonly)
       open (itube,file=tubename,status='old',action='read')
       call getenv_loc('IONPAR',ionname)
       open (iion,file=ionname,status='old',action='read')
+      call getenv_loc('IONPAR_NUCL',ionnuclname)
+      open (iionnucl,file=ionnuclname,status='old',action='read')
 
 #ifndef CRYST_SC
       call getenv_loc('ROTPARPDB',rotname_pdb)