adding ebend_nucl to UCGM+some further reading
[unres4.git] / source / unres / io_config.f90
index 5153cc9..25a5e24 100644 (file)
 !
       allocate(dsc(ntyp1)) !(ntyp1)
       allocate(dsc_inv(ntyp1)) !(ntyp1)
+      allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
+      allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
+      allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
       allocate(nbondterm(ntyp)) !(ntyp)
       allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
       allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
-      allocate(msc(ntyp+1)) !(ntyp+1)
-      allocate(isc(ntyp+1)) !(ntyp+1)
-      allocate(restok(ntyp+1)) !(ntyp+1)
+      allocate(msc(ntyp+1,5)) !(ntyp+1)
+      allocate(isc(ntyp+1,5)) !(ntyp+1)
+      allocate(restok(ntyp+1,5)) !(ntyp+1)
       allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
       allocate(long_r_sidechain(ntyp))
       allocate(short_r_sidechain(ntyp))
         endif
       enddo
 #else
-      read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp,ip,pstok
-      do i=1,ntyp
+      read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
+      do i=1,ntyp_molec(1)
         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
-         j=1,nbondterm(i)),msc(i),isc(i),restok(i)
+         j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
         dsc(i) = vbldsc0(1,i)
         if (i.eq.10) then
           dsc_inv(i)=0.0D0
         endif
       enddo
 #endif
+      read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
+      do i=1,ntyp_molec(2)
+        nbondterm_nucl(i)=1
+        read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
+!        dsc(i) = vbldsc0_nucl(1,i)
+!        if (i.eq.10) then
+!          dsc_inv(i)=0.0D0
+!        else
+!          dsc_inv(i)=1.0D0/dsc(i)
+!        endif
+      enddo
+!      read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
+!      do i=1,ntyp_molec(2)
+!        read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),& 
+!        aksc_nucl(j,i),abond0_nucl(j,i),&
+!         j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
+!        dsc(i) = vbldsc0(1,i)
+!        if (i.eq.10) then
+!          dsc_inv(i)=0.0D0
+!        else
+!          dsc_inv(i)=1.0D0/dsc(i)
+!        endif
+!      enddo
+
       if (lprint) then
         write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
          'inertia','Pstok'
-        write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
+        write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
         do i=1,ntyp
           write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
-            vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
+            vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
           do j=2,nbondterm(i)
             write (iout,'(13x,3f10.5)') &
               vbldsc0(j,i),aksc(j,i),abond0(j,i)
       close (ithep_pdb)
 #endif
       close(ithep)
+!--------------- Reading theta parameters for nucleic acid-------
+      read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
+      ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
+      nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
+      allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
+      allocate(aa0thet_nucl(nthetyp_nucl+1,&
+        nthetyp_nucl+1,nthetyp_nucl+1))
+!(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
+      allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
+        nthetyp_nucl+1,nthetyp_nucl+1))
+!(maxtheterm,-maxthetyp1:maxthetyp1,&
+!        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
+      allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+        nthetyp_nucl+1,nthetyp_nucl+1))
+      allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+        nthetyp_nucl+1,nthetyp_nucl+1))
+      allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+        nthetyp_nucl+1,nthetyp_nucl+1))
+      allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+        nthetyp_nucl+1,nthetyp_nucl+1))
+!(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
+!        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
+      allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
+         nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
+      allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
+         nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
+
+!(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
+!        -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+
+      read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
+
+      aa0thet_nucl(:,:,:)=0.0d0
+      aathet_nucl(:,:,:,:)=0.0d0
+      bbthet_nucl(:,:,:,:,:)=0.0d0
+      ccthet_nucl(:,:,:,:,:)=0.0d0
+      ddthet_nucl(:,:,:,:,:)=0.0d0
+      eethet_nucl(:,:,:,:,:)=0.0d0
+      ffthet_nucl(:,:,:,:,:,:)=0.0d0
+      ggthet_nucl(:,:,:,:,:,:)=0.0d0
+
+      do i=1,nthetyp_nucl
+        do j=1,nthetyp_nucl
+          do k=1,nthetyp_nucl
+            read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
+            read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
+            read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
+            read (ithep_nucl,*,end=111,err=111) &
+            (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
+            (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
+            (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
+            (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
+            read (ithep_nucl,*,end=111,err=111) &
+           (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
+              ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
+              llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
+          enddo
+        enddo
+      enddo
 
 !-------------------------------------------
       allocate(nlob(ntyp1)) !(ntyp1)
       character(len=80) :: card
       real(kind=8),dimension(3,20) :: sccor
       integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin      !rescode,
-      integer :: isugar
+      integer :: isugar,molecprev
       character*1 :: sugar
       real(kind=8) :: cou
       real(kind=8),dimension(3) :: ccc
 !              nres_molec(molecule)=nres_molec(molecule)+1
             else
              molecule=2
-              itype(ires,molecule)=rescode(ires,res(3:4),0,molecule)
+              itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
 !              nres_molec(molecule)=nres_molec(molecule)+1
             endif
             endif
                    atom.ne.'OXT' .and. atom(:2).ne.'3H' &
                    .and. atom.ne.'P  '.and. &
                   atom(1:1).ne.'H' .and. &
-                  atom.ne.'OP1' .and. atom.ne.'OP2 ') then
+                  atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
+                  ) then
 !            write (iout,*) "sidechain ",atom
 !            write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
                  if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
             read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
               endif
           endif
-        endif
+        else if ((ions).and.(card(1:6).eq.'HETATM')) then
+         
+          read (card(12:16),*) atom
+          print *,"HETATOM", atom
+          read (card(18:20),'(a3)') res
+          if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
+          (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG')           &
+          .or.(atom(1:2).eq.'K ')) &
+          then
+           ires=ires+1
+           if (molecule.ne.5) molecprev=molecule
+           molecule=5
+           nres_molec(molecule)=nres_molec(molecule)+1
+           itype(ires,molecule)=rescode(ires,res(2:4),0,molecule)
+           read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+          endif
+        endif !atom
       enddo
    10 write (iout,'(a,i5)') ' Number of residues found: ',ires
       if (ires.eq.0) return
       if (ires_old.ne.ires) nres_molec(molecule)=nres_molec(molecule)-2
 !      print *,'I have', nres_molec(:)
       
-      do k=1,5 
+      do k=1,4 ! ions are without dummy 
        if (nres_molec(k).eq.0) cycle
       do i=2,nres-1
 !        write (iout,*) i,itype(i,1)
       nsup=nres
       nstart_sup=1
 !      print *,"molecule",molecule
-      if (itype(nres,1).ne.10) then
+      if ((itype(nres,1).ne.10)) then
         nres=nres+1
+          if (molecule.eq.5) molecule=molecprev
         itype(nres,molecule)=ntyp1_molec(molecule)
         nres_molec(molecule)=nres_molec(molecule)+1
         if (unres_pdb) then
 ! NOW LETS ROCK! SORTING
       allocate(c_temporary(3,2*nres))
       allocate(itype_temporary(nres,5))
+      allocate(molnum(nres))
        itype_temporary(:,:)=0
       seqalingbegin=1
       do k=1,5
 
           enddo
           itype_temporary(seqalingbegin,k)=itype(i,k)
+          molnum(seqalingbegin)=k
           seqalingbegin=seqalingbegin+1
          endif
         enddo
         if(me.eq.king.or..not.out1file) &
          write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
           eta
-        gamp=scal_fric*(pstok+rwat)*eta
+        gamp=scal_fric*(pstok(i)+rwat)*eta
         stdfp=dsqrt(2*Rb*t_bath/d_time)
         allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
         do i=1,ntyp
-          gamsc(i)=scal_fric*(restok(i)+rwat)*eta  
+          gamsc(i)=scal_fric*(restok(i,1)+rwat)*eta  
           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
         enddo 
         if(me.eq.king.or..not.out1file)then
          " stochastic forces of fully exposed sites"
          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
          do i=1,ntyp
-          write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i),&
+          write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
            gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
          enddo
         endif
 ! Get parameter filenames and open the parameter files.
       call getenv_loc('BONDPAR',bondname)
       open (ibond,file=bondname,status='old',readonly,shared)
+      call getenv_loc('BONDPAR_NUCL',bondname_nucl)
+      open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
       call getenv_loc('THETPAR',thetname)
       open (ithep,file=thetname,status='old',readonly,shared)
       call getenv_loc('ROTPAR',rotname)
       open (ielep,file=elename,status='old',readonly,shared)
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old',readonly,shared)
+
+      call getenv_loc('THETPAR_NUCL',thetname_nucl)
+      open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
+      call getenv_loc('ROTPAR_NUCL',rotname_nucl)
+      open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
+      call getenv_loc('TORPAR_NUCL',torname_nucl)
+      open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
+      call getenv_loc('TORDPAR_NUCL',tordname_nucl)
+      open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
+      call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
+      open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
+
+
 #elif (defined CRAY) || (defined AIX)
       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
         action='read')
 ! Get parameter filenames and open the parameter files.
       call getenv_loc('BONDPAR',bondname)
       open (ibond,file=bondname,status='old',action='read')
+      call getenv_loc('BONDPAR_NUCL',bondname_nucl)
+      open (ibond_nucl,file=bondname_nucl,status='old',action='read')
+
 !      print *,"Processor",myrank," opened file IBOND" 
       call getenv_loc('THETPAR',thetname)
       open (ithep,file=thetname,status='old',action='read')
 !      print *,"Processor",myrank," opened file IELEP" 
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old',action='read')
+
+      call getenv_loc('THETPAR_NUCL',thetname_nucl)
+      open (ithep_nucl,file=thetname_nucl,status='old',action='read')
+      call getenv_loc('ROTPAR_NUCL',rotname_nucl)
+      open (irotam_nucl,file=rotname_nucl,status='old',action='read')
+      call getenv_loc('TORPAR_NUCL',torname_nucl)
+      open (itorp_nucl,file=torname_nucl,status='old',action='read')
+      call getenv_loc('TORDPAR_NUCL',tordname_nucl)
+      open (itordp_nucl,file=tordname_nucl,status='old',action='read')
+      call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
+      open (isidep_nucl,file=sidename_nucl,status='old',action='read')
+
       call getenv_loc('LIPTRANPAR',liptranname)
       open (iliptranpar,file=liptranname,status='old',action='read')
       call getenv_loc('TUBEPAR',tubename)
 ! Get parameter filenames and open the parameter files.
       call getenv_loc('BONDPAR',bondname)
       open (ibond,file=bondname,status='old')
+      call getenv_loc('BONDPAR_NUCL',bondname_nucl)
+      open (ibond_nucl,file=bondname_nucl,status='old')
+
       call getenv_loc('THETPAR',thetname)
       open (ithep,file=thetname,status='old')
       call getenv_loc('ROTPAR',rotname)
       open (ielep,file=elename,status='old')
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old')
+
+      open (ithep_nucl,file=thetname_nucl,status='old')
+      call getenv_loc('ROTPAR_NUCL',rotname_nucl)
+      open (irotam_nucl,file=rotname_nucl,status='old')
+      call getenv_loc('TORPAR_NUCL',torname_nucl)
+      open (itorp_nucl,file=torname_nucl,status='old')
+      call getenv_loc('TORDPAR_NUCL',tordname_nucl)
+      open (itordp_nucl,file=tordname_nucl,status='old')
+      call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
+      open (isidep_nucl,file=sidename_nucl,status='old')
+
       call getenv_loc('LIPTRANPAR',liptranname)
       open (iliptranpar,file=liptranname,status='old')
       call getenv_loc('TUBEPAR',tubename)
 ! Get parameter filenames and open the parameter files.
       call getenv_loc('BONDPAR',bondname)
       open (ibond,file=bondname,status='old',action='read')
+      call getenv_loc('BONDPAR_NUCL',bondname_nucl)
+      open (ibond_nucl,file=bondname_nucl,status='old',action='read')
       call getenv_loc('THETPAR',thetname)
       open (ithep,file=thetname,status='old',action='read')
       call getenv_loc('ROTPAR',rotname)
       open (ielep,file=elename,status='old',readonly)
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old',readonly)
+
+      call getenv_loc('THETPAR_NUCL',thetname_nucl)
+      open (ithep_nucl,file=thetname_nucl,status='old',action='read')
+      call getenv_loc('ROTPAR_NUCL',rotname_nucl)
+      open (irotam_nucl,file=rotname_nucl,status='old',action='read')
+      call getenv_loc('TORPAR_NUCL',torname_nucl)
+      open (itorp_nucl,file=torname_nucl,status='old',action='read')
+      call getenv_loc('TORDPAR_NUCL',tordname_nucl)
+      open (itordp_nucl,file=tordname_nucl,status='old',action='read')
+      call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
+      open (isidep_nucl,file=sidename_nucl,status='old',action='read')
+
       call getenv_loc('LIPTRANPAR',liptranname)
       open (iliptranpar,file=liptranname,status='old',action='read')
       call getenv_loc('TUBEPAR',tubename)