working microcanonical for CA2+ protein

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

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)