working microcanonical for CA2+ protein
[unres4.git] / source / unres / io_config.f90
index 6ae8cfd..8e913d0 100644 (file)
       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,&
           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)
       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
              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
             do j=1,3
               c(j,ires)=c(j,ires)+ccc(j)/5.0
             enddo
+             print *,counter,molecule
              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
               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
            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
 ! 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 
 
           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
 !      character(len=80) :: ucase
       character(len=320) :: controlcard
 !el local variables
-      integer :: i
+      integer :: i,j
       real(kind=8) :: eta
 
       call card_concat(controlcard,.true.)
         if(me.eq.king.or..not.out1file) &
          write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
           eta
-        gamp=scal_fric*(pstok(i)+rwat)*eta
-        stdfp=dsqrt(2*Rb*t_bath/d_time)
-        allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
+!        allocate(gamp
+        do j=1,5 !types of molecules
+        gamp(j)=scal_fric*(pstok(j)+rwat)*eta
+        stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
+        enddo
+        allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
+        do j=1,5 !types of molecules
         do i=1,ntyp
-          gamsc(i)=scal_fric*(restok(i,1)+rwat)*eta  
-          stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
+          gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta  
+          stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
         enddo 
+        enddo
+
         if(me.eq.king.or..not.out1file)then
          write (iout,'(/2a/)') &
          "Radii of site types and friction coefficients and std's of",&
          " stochastic forces of fully exposed sites"
-         write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
+         write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
          do i=1,ntyp
           write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
-           gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
+           gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
          enddo
         endif
       else if (tbf) then
       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" 
       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)
       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)