working martini
[unres4.git] / source / unres / io_config.F90
index 9b036aa..cdd6231 100644 (file)
           dsc_inv(i)=1.0D0/dsc(i)
         endif
       enddo
+!      ip(1)=0.0001d0
+!      isc(:,1)=0.0001d0
 #endif
       read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
       do i=1,ntyp_molec(2)
              read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
              print *,msc(i,5),restok(i,5)
             enddo
-            ip(5)=0.2
+!            ip(5)=0.2
 !            isc(5)=0.2
             read (iion,*) ncatprotparm
             allocate(catprm(ncatprotparm,4))
          enddo  
        endif
       enddo
+#ifdef SC_END
+      allocate(nterm_scend(2,ntyp))
+      allocate(arotam_end(0:6,2,ntyp))
+      nterm_scend=0
+      arotam_end=0.0d0
+      read (irotam_end,*) ijunk
+!c      write (iout,*) "ijunk",ijunk
+      do i=1,ntyp
+        if (i.eq.10) cycle
+        do j=1,2
+          read (irotam_end,'(a)')
+          read (irotam_end,*) nterm_scend(j,i)
+!c          write (iout,*) "i",i," j",j," nterm",nterm_scend(j,i)
+          do k=0,nterm_scend(j,i)
+            read (irotam_end,*) ijunk,arotam_end(k,j,i)
+!c            write (iout,*) "k",k," arotam",arotam_end(k,j,i)
+          enddo
+        enddo
+      enddo
+!c      lprint=.true.
+      if (lprint) then
+        write (iout,'(a)') &
+         "Parameters of the local potentials of sidechain ends"
+        do i=1,ntyp
+          write (iout,'(5x,9x,2hp-,a3,6x,9x,a3,2h-p)')&
+          restyp(i,1),restyp(i,1)
+          do j=0,max0(nterm_scend(1,i),nterm_scend(2,i))
+            write (iout,'(i5,2f20.10)') &
+             j,arotam_end(j,1,i),arotam_end(j,2,i)
+          enddo
+        enddo
+      endif
+!c      lprint=.false.
+#endif
+
 !---------reading nucleic acid parameters for rotamers-------------------
       allocate(sc_parmin_nucl(9,ntyp_molec(2)))      !(maxsccoef,ntyp)
       do i=1,ntyp_molec(2)
         sigiso2(i,j)=sigiso2(j,i)
 !        print *,"ATU",sigma(j,i),sigma(i,j),i,j
         nstate(i,j) = nstate(j,i)
-        dtail(1,i,j) = dtail(1,j,i)
-        dtail(2,i,j) = dtail(2,j,i)
+        dtail(1,i,j) = dtail(2,j,i)
+        dtail(2,i,j) = dtail(1,j,i)
         DO k = 1, 4
          alphasur(k,i,j) = alphasur(k,j,i)
          wstate(k,i,j) = wstate(k,j,i)
       !HERE THE MASS of MARTINI
       write(*,*) "before MARTINI PARAM"
       do i=1,ntyp_molec(4)
-       msc(i,4)=0.0d0
-       mp(4)=72.0d0
+       msc(i,4)=72.0d0
+       mp(4)=0.0d0
        isc(i,4)=0.d0
       enddo
       ip(4)=0.0
+      msc(ntyp_molec(4)+1,4)=0.1d0
       !relative dielectric constant = 15 for implicit screening
       k_coulomb_lip=332.0d0/15.0d0
       !kbond = 1250 kJ/(mol*nm*2)
 
 ! Ions by Aga
 
-       allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
-       allocate(alphapolcat2(ntyp,ntyp))
-       allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
-       allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
-       allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(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))
+       allocate(alphapolcat(ntyp,-1:ntyp_molec(5)),epsheadcat(ntyp,-1:ntyp_molec(5)),sig0headcat(ntyp,-1:ntyp_molec(5)))
+       allocate(alphapolcat2(ntyp,-1:ntyp_molec(5)))
+       allocate(sigiso1cat(ntyp,-1:ntyp_molec(5)),rborn1cat(ntyp,-1:ntyp_molec(5)),rborn2cat(ntyp,-1:ntyp_molec(5)),sigmap1cat(ntyp,-1:ntyp_molec(5)))
+       allocate(sigmap2cat(ntyp,-1:ntyp_molec(5)),sigiso2cat(ntyp,-1:ntyp_molec(5)))
+       allocate(chis1cat(ntyp,-1:ntyp_molec(5)),chis2cat(ntyp,-1:ntyp_molec(5)),wquadcat(ntyp,-1:ntyp_molec(5)),chipp1cat(ntyp,-1:ntyp_molec(5)),chipp2cat(ntyp,-1:ntyp_molec(5)))
+       allocate(epsintabcat(ntyp,-1:ntyp_molec(5)))
+       allocate(dtailcat(2,ntyp,-1:ntyp_molec(5)))
+       allocate(alphasurcat(4,ntyp,-1:ntyp_molec(5)),alphisocat(4,ntyp,-1:ntyp_molec(5)))
+       allocate(wqdipcat(2,ntyp,-1:ntyp_molec(5)))
+       allocate(wstatecat(4,ntyp,-1:ntyp_molec(5)))
+       allocate(dheadcat(2,2,ntyp,-1:ntyp_molec(5)))
+       allocate(nstatecat(ntyp,-1:ntyp_molec(5)))
+       allocate(debaykapcat(ntyp,-1:ntyp_molec(5)))
+
+      if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,-1:ntyp1))
+      if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,-1:ntyp1))
 !      if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
-      if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
-      if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
+      if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,-1:ntyp1)) !(ntyp,ntyp)
+      if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,-1:ntyp1)) !(ntyp,ntyp)
 
 
             if (.not.allocated(ichargecat))&
              read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
              print *,msc(i,5),restok(i,5)
             enddo
-            ip(5)=0.2
+           ! ip(5)=0.2
            ! mp(5)=0.2
              pstok(5)=3.0
 !DIR$ NOUNROLL 
-      do j=1,ntyp_molec(5)-1 ! this is without Zn will be modified for ALL tranistion metals
+      do j=-1,ntyp_molec(5)-1 ! this is without Zn will be modified for ALL tranistion metals
+       if (j.eq.0) cycle
        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), &
-       chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
+       chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), & !6
 
-       (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
+       (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&!12
 !       chiscat(i,j),chiscat(j,i), &
        chis1cat(i,j),chis2cat(i,j), &
 
-       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),&
+       nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !19                          !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),&!23
        dtailcat(1,i,j),dtailcat(2,i,j), &
-       epsheadcat(i,j),sig0headcat(i,j), &
+       epsheadcat(i,j),sig0headcat(i,j), &!27
 !wdipcat = w1 , w2
 !       rborncat(i,j),rborncat(j,i),&
        rborn1cat(i,j),rborn2cat(i,j),&
-       (wqdipcat(k,i,j),k=1,2), &
-       alphapolcat(i,j),alphapolcat2(j,i), &
+       (wqdipcat(k,i,j),k=1,2), &!31
+       alphapolcat(i,j),alphapolcat2(j,i), &!33
        (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
 
        if (chi1cat(i,j).gt.0.9) write (*,*) "WTF ANISO", i,j, chi1cat(i,j)
 
        END DO
       END DO
-      allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
+      allocate(aa_aq_cat(-ntyp:ntyp,-1:ntyp_molec(5)),&
+               bb_aq_cat(-ntyp:ntyp,-1:ntyp_molec(5)))
       do i=1,ntyp
-        do j=1,ntyp_molec(5)
+        do j=-1,ntyp_molec(5)
+          if (j.eq.0) cycle
           epsij=epscat(i,j)
           rrij=sigmacat(i,j)
           rrij=rrij**expon
       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
         ' v3ss:',v3ss
       endif
+
+!------------MARTINI-PROTEIN-parameters-------------------------
+       allocate(alphapolmart(ntyp,ntyp),epsheadmart(ntyp,ntyp_molec(4)),sig0headmart(ntyp,ntyp_molec(4)))
+       allocate(alphapolmart2(ntyp,ntyp))
+       allocate(sigiso1mart(ntyp,ntyp_molec(4)),rborn1mart(ntyp,ntyp_molec(4)),rborn2mart(ntyp,ntyp_molec(4)),sigmap1mart(ntyp,ntyp_molec(4)))
+       allocate(sigmap2mart(ntyp,ntyp_molec(4)),sigiso2mart(ntyp,ntyp_molec(4)))
+       allocate(chis1mart(ntyp,ntyp_molec(4)),chis2mart(ntyp,ntyp_molec(4)),wquadmart(ntyp,ntyp_molec(4)),chipp1mart(ntyp,ntyp_molec(4)),chipp2mart(ntyp,ntyp_molec(4)))
+       allocate(epsintabmart(ntyp,ntyp_molec(4)))
+       allocate(dtailmart(2,ntyp,ntyp_molec(4)))
+       allocate(alphasurmart(4,ntyp,ntyp_molec(4)),alphisomart(4,ntyp,ntyp_molec(4)))
+       allocate(wqdipmart(2,ntyp,ntyp_molec(4)))
+       allocate(wstatemart(4,ntyp,ntyp_molec(4)))
+       allocate(dheadmart(2,2,ntyp,ntyp_molec(4))) 
+       allocate(nstatemart(ntyp,ntyp_molec(4)))
+       allocate(debaykapmart(ntyp,ntyp_molec(4)))
+
+      if (.not.allocated(epsmart)) allocate (epsmart(0:ntyp1,ntyp1))
+      if (.not.allocated(sigmamart)) allocate(sigmamart(0:ntyp1,ntyp1))
+!      if (.not.allocated(chimart)) allomarte(chimart(ntyp1,ntyp1)) !(ntyp,ntyp)
+      if (.not.allocated(chi1mart)) allocate(chi1mart(ntyp1,ntyp1)) !(ntyp,ntyp)
+      if (.not.allocated(chi2mart)) allocate(chi2mart(ntyp1,ntyp1)) !(ntyp,ntyp)
+
+!DIR$ NOUNROLL 
+      do i=1,ntyp-3 ! there are phosporylated missing
+      do j=1,ntyp_molec(4) ! this is without Zn will be modified for ALL tranistion metals
+!       do j=1,ntyp_molec(5) 
+        print *,"lipmart",i,j
+!        write (*,*) "Im in ALAB", i, " ", j
+        read(imartprot,*) &
+       epsmart(i,j),sigmamart(i,j), &
+!       chimart(i,j),chimart(j,i),chippmart(i,j),chippmart(j,i), &
+       chi1mart(i,j),chi2mart(i,j),chipp1mart(i,j),chipp2mart(i,j), & !6
+
+       (alphasurmart(k,i,j),k=1,4),sigmap1mart(i,j),sigmap2mart(i,j),&!12
+!       chismart(i,j),chismart(j,i), &
+       chis1mart(i,j),chis2mart(i,j), &
+
+       nstatemart(i,j),(wstatemart(k,i,j),k=1,4), & !19                          !5 w tej lini - 1 integer pierwszy
+       dheadmart(1,1,i,j),dheadmart(1,2,i,j),dheadmart(2,1,i,j),dheadmart(2,2,i,j),&!23
+       dtailmart(1,i,j),dtailmart(2,i,j), &
+       epsheadmart(i,j),sig0headmart(i,j), &!27 
+!wdipmart = w1 , w2
+!       rbornmart(i,j),rbornmart(j,i),&
+       rborn1mart(i,j),rborn2mart(i,j),&
+       (wqdipmart(k,i,j),k=1,2), &!31
+       alphapolmart(i,j),alphapolmart2(j,i), &!33
+       (alphisomart(k,i,j),k=1,4),sigiso1mart(i,j),sigiso2mart(i,j),epsintabmart(i,j),debaykapmart(i,j) 
+       enddo
+       enddo
+      allocate(aa_aq_mart(-ntyp:ntyp,ntyp_molec(4)),&
+               bb_aq_mart(-ntyp:ntyp,ntyp_molec(4)))
+      do i=1,ntyp-3 ! still no phophorylated residues
+        do j=1,ntyp_molec(4)
+          if (j.eq.0) cycle
+          epsij=epsmart(i,j)
+          rrij=sigmamart(i,j)
+          rrij=rrij**expon
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_aq_mart(i,j)=epsij*rrij*rrij
+          bb_aq_mart(i,j)=-sigeps*epsij*rrij
+         enddo
+       enddo
+
+
+
+
+
+
+
+
+
+
+
       if (shield_mode.gt.0) then
       pi=4.0D0*datan(1.0D0)
 !C VSolvSphere the volume of solving sphere
       cou=1
         write (iout,*) "symetr", symetr
       do i=1,nres
-      lll=lll+1
+       lll=lll+1
 !      write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
-      if (i.gt.1) then
-      if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
-      chain_length=lll-1
-      kkk=kkk+1
+!      if (i.gt.1) then
+!      if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
+!      chain_length=lll-1
+!      kkk=kkk+1
 !       write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
-      lll=1
-      endif
-      endif
+!      lll=1
+!      endif
+!      endif
         do j=1,3
           cref(j,i,cou)=c(j,i)
           cref(j,i+nres,cou)=c(j,i+nres)
           endif
          enddo
       enddo
-      write (iout,*) chain_length
-      if (chain_length.eq.0) chain_length=nres
-      do j=1,3
-      chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
-      chain_rep(j,chain_length+nres,symetr) &
-      =chain_rep(j,chain_length+nres,1)
-      enddo
 ! diagnostic
 !       write (iout,*) "spraw lancuchy",chain_length,symetr
 !       do i=1,4
       dc(j,0)=c(j,1)
       enddo
 
-      if (symetr.gt.1) then
-       call permut(symetr)
-       nperm=1
-       do i=1,symetr
-       nperm=nperm*i
-       enddo
-       do i=1,nperm
-       write(iout,*) (tabperm(i,kkk),kkk=1,4)
-       enddo
-       do i=1,nperm
-        cou=0
-        do kkk=1,symetr
-         icha=tabperm(i,kkk)
-         write (iout,*) i,icha
-         do lll=1,chain_length
-          cou=cou+1
-           if (cou.le.nres) then
-           do j=1,3
-            kupa=mod(lll,chain_length)
-            iprzes=(kkk-1)*chain_length+lll
-            if (kupa.eq.0) kupa=chain_length
-            write (iout,*) "kupa", kupa
-            cref(j,iprzes,i)=chain_rep(j,kupa,icha)
-            cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
-          enddo
-          endif
-         enddo
-        enddo
-       enddo
-       endif
+!      if (symetr.gt.1) then
+!       call permut(symetr)
+!       nperm=1
+!       do i=1,symetr
+!       nperm=nperm*i
+!       enddo
+!      do i=1,nperm
+!      write(iout,*) (tabperm(i,kkk),kkk=1,4)
+!      enddo
+!      do i=1,nperm
+!       cou=0
+!       do kkk=1,symetr
+!        icha=tabperm(i,kkk)
+!        write (iout,*) i,icha
+!        do lll=1,chain_length
+!         cou=cou+1
+!          if (cou.le.nres) then
+!          do j=1,3
+!           kupa=mod(lll,chain_length)
+!           iprzes=(kkk-1)*chain_length+lll
+!           if (kupa.eq.0) kupa=chain_length
+!           write (iout,*) "kupa", kupa
+!           cref(j,iprzes,i)=chain_rep(j,kupa,icha)
+!           cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
+!         enddo
+!         endif
+!        enddo
+!       enddo
+!      enddo
+!      endif
 !-koniec robienia kopii
 ! diag
       do kkk=1,nperm
 ! CUTOFFF ON ELECTROSTATICS
       call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
-      write(iout,*) "R_CUT_ELE=",r_cut_ele
+      
+      write(iout,*) "R_CUT_ELE=",r_cut_ele,rlamb_ele
       call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0)
       call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0)
       call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0)
 
+      call reada(controlcard,"DELTA",graddelta,1.0d-5)
 ! Lipidec parameters
       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
 !      print *,"Processor",myrank," opened file ITHEP" 
       call getenv_loc('ROTPAR',rotname)
       open (irotam,file=rotname,status='old',action='read')
+#ifdef SC_END
+      call getenv_loc('ROTPAR_END',rotname_end)
+      open (irotam_end,file=rotname_end,status='old',action='read')
+#endif
 !      print *,"Processor",myrank," opened file IROTAM" 
       call getenv_loc('TORPAR',torname)
       open (itorp,file=torname,status='old',action='read')
       open (ithep,file=thetname,status='old')
       call getenv_loc('ROTPAR',rotname)
       open (irotam,file=rotname,status='old')
+#ifdef SC_END
+      call getenv_loc('ROTPAR_END',rotname_end)
+      open (irotam_end,file=rotname_end,status='old')
+#endif
       call getenv_loc('TORPAR',torname)
       open (itorp,file=torname,status='old')
       call getenv_loc('TORDPAR',tordname)
       open (ithep,file=thetname,status='old',action='read')
       call getenv_loc('ROTPAR',rotname)
       open (irotam,file=rotname,status='old',action='read')
+#ifdef SC_END
+      call getenv_loc('ROTPAR_END',rotname_end)
+      open (irotam_end,file=rotname_end,status='old',action='read')
+#endif
       call getenv_loc('TORPAR',torname)
       open (itorp,file=torname,status='old',action='read')
       call getenv_loc('TORDPAR',tordname)
       open (ilipbond,file=lipbondname,status='old',action='read')
       call getenv_loc('LIPNONBOND',lipnonbondname)
       open (ilipnonbond,file=lipnonbondname,status='old',action='read')
+      call getenv_loc('LIPPROT',lipprotname)
+      open (imartprot,file=lipprotname,status='old',action='read')
+
       call getenv_loc('TUBEPAR',tubename)
       open (itube,file=tubename,status='old',action='read')
       call getenv_loc('IONPAR',ionname)