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 a66638b..2906fb5 100644 (file)
--- a/source/unres/io_config.F90
+++ b/source/unres/io_config.F90
@@ -3250,14 +3250,21 @@
        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)
+             read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
              print *,msc(i,5),restok(i,5)
             enddo
             ip(5)=0.2
@@ -3280,6 +3287,17 @@
 !       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
 
       
@@ -4294,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'))