working gradient for cations

[unres4.git] / source / unres / io_config.F90

diff --git a/source/unres/io_config.F90 b/source/unres/io_config.F90
index 0a1e048..2906fb5 100644 (file)
--- a/source/unres/io_config.F90
+++ b/source/unres/io_config.F90
@@ -867,6 +867,7 @@
           enddo
         enddo
       endif
+       if (oldion.eq.1) then
             do i=1,ntyp_molec(5)
              read(iion,*) msc(i,5),restok(i,5)
              print *,msc(i,5),restok(i,5)
@@ -879,6 +880,7 @@
             read (iion,*)  (catprm(i,k),i=1,ncatprotparm)
             enddo
             print *, catprm
+         endif
 !            read (iion,*) (vcatprm(k),k=1,ncatprotpram)
 !----------------------------------------------------
       allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
@@ -2730,6 +2732,8 @@
        END DO
       END DO
       endif
+
+
 !      goto 50
 !--------------------- GBV potential -----------------------------------
        case(5)
@@ -3230,6 +3234,72 @@
 !      v1ss=0.0d0
 !      v2ss=0.0d0
 !      v3ss=0.0d0
+
+! Ions by Aga
+
+       allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
+       allocate(sigiso1cat(ntyp,ntyp),rborncat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
+       allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
+       allocate(chiscat(ntyp,ntyp),wquadcat(ntyp,ntyp),chippcat(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)
+      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
+
+      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), &
+       (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
+       chiscat(i,j),chiscat(j,i), &
+       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), &
+       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)
+          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 
+      endif
+
       
       if(me.eq.king) then
       write (iout,'(/a)') "Disulfide bridge parameters:"
@@ -4242,6 +4312,7 @@
       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'))
@@ -4259,6 +4330,8 @@
 ! elemode = 0 is orignal UNRES electrostatics
 ! elemode = 1 is "Momo" potentials in progress
 ! elemode = 2 is in development EVALD
+
+
       write (iout,*) TUBEmode,"TUBEMODE"
       if (TUBEmode.gt.0) then
        call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)