critical bug in set_matrices
[unres4.git] / source / unres / io_config.f90
index 8e913d0..378086b 100644 (file)
       character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
       logical :: lprint,LaTeX
       real(kind=8),dimension(3,3,maxlob) :: blower     !(3,3,maxlob)
-      real(kind=8),dimension(13) :: b
+      real(kind=8),dimension(13) :: buse
       character(len=3) :: lancuch      !,ucase
 !el  local variables
       integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
 
 #endif      
       if (lprint) then
+        l=3
         write (iout,'(/a/)') 'Torsional constants:'
         do i=1,nsccortyp
           do j=1,nsccortyp
 
       do i=0,nloctyp-1
         read (ifourier,*,end=115,err=115)
-        read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
+        read (ifourier,*,end=115,err=115) (buse(ii),ii=1,13)
         if (lprint) then
           write (iout,*) 'Type',i
-          write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
+          write (iout,'(a,i2,a,f10.5)') ('buse(',ii,')=',buse(ii),ii=1,13)
         endif
-        B1(1,i)  = b(3)
-        B1(2,i)  = b(5)
-        B1(1,-i) = b(3)
-        B1(2,-i) = -b(5)
-!        b1(1,i)=0.0d0
-!        b1(2,i)=0.0d0
-        B1tilde(1,i) = b(3)
-        B1tilde(2,i) =-b(5)
-        B1tilde(1,-i) =-b(3)
-        B1tilde(2,-i) =b(5)
-!        b1tilde(1,i)=0.0d0
-!        b1tilde(2,i)=0.0d0
-        B2(1,i)  = b(2)
-        B2(2,i)  = b(4)
-        B2(1,-i)  =b(2)
-        B2(2,-i)  =-b(4)
-
-!        b2(1,i)=0.0d0
-!        b2(2,i)=0.0d0
-        CC(1,1,i)= b(7)
-        CC(2,2,i)=-b(7)
-        CC(2,1,i)= b(9)
-        CC(1,2,i)= b(9)
-        CC(1,1,-i)= b(7)
-        CC(2,2,-i)=-b(7)
-        CC(2,1,-i)=-b(9)
-        CC(1,2,-i)=-b(9)
+        B1(1,i)  = buse(3)
+        B1(2,i)  = buse(5)
+        B1(1,-i) = buse(3)
+        B1(2,-i) = -buse(5)
+!        buse1(1,i)=0.0d0
+!        buse1(2,i)=0.0d0
+        B1tilde(1,i) = buse(3)
+        B1tilde(2,i) =-buse(5)
+        B1tilde(1,-i) =-buse(3)
+        B1tilde(2,-i) =buse(5)
+!        buse1tilde(1,i)=0.0d0
+!        buse1tilde(2,i)=0.0d0
+        B2(1,i)  = buse(2)
+        B2(2,i)  = buse(4)
+        B2(1,-i)  =buse(2)
+        B2(2,-i)  =-buse(4)
+
+!        buse2(1,i)=0.0d0
+!        buse2(2,i)=0.0d0
+        CC(1,1,i)= buse(7)
+        CC(2,2,i)=-buse(7)
+        CC(2,1,i)= buse(9)
+        CC(1,2,i)= buse(9)
+        CC(1,1,-i)= buse(7)
+        CC(2,2,-i)=-buse(7)
+        CC(2,1,-i)=-buse(9)
+        CC(1,2,-i)=-buse(9)
 !        CC(1,1,i)=0.0d0
 !        CC(2,2,i)=0.0d0
 !        CC(2,1,i)=0.0d0
 !        CC(1,2,i)=0.0d0
-        Ctilde(1,1,i)=b(7)
-        Ctilde(1,2,i)=b(9)
-        Ctilde(2,1,i)=-b(9)
-        Ctilde(2,2,i)=b(7)
-        Ctilde(1,1,-i)=b(7)
-        Ctilde(1,2,-i)=-b(9)
-        Ctilde(2,1,-i)=b(9)
-        Ctilde(2,2,-i)=b(7)
+        Ctilde(1,1,i)=buse(7)
+        Ctilde(1,2,i)=buse(9)
+        Ctilde(2,1,i)=-buse(9)
+        Ctilde(2,2,i)=buse(7)
+        Ctilde(1,1,-i)=buse(7)
+        Ctilde(1,2,-i)=-buse(9)
+        Ctilde(2,1,-i)=buse(9)
+        Ctilde(2,2,-i)=buse(7)
 
 !        Ctilde(1,1,i)=0.0d0
 !        Ctilde(1,2,i)=0.0d0
 !        Ctilde(2,1,i)=0.0d0
 !        Ctilde(2,2,i)=0.0d0
-        DD(1,1,i)= b(6)
-        DD(2,2,i)=-b(6)
-        DD(2,1,i)= b(8)
-        DD(1,2,i)= b(8)
-        DD(1,1,-i)= b(6)
-        DD(2,2,-i)=-b(6)
-        DD(2,1,-i)=-b(8)
-        DD(1,2,-i)=-b(8)
+        DD(1,1,i)= buse(6)
+        DD(2,2,i)=-buse(6)
+        DD(2,1,i)= buse(8)
+        DD(1,2,i)= buse(8)
+        DD(1,1,-i)= buse(6)
+        DD(2,2,-i)=-buse(6)
+        DD(2,1,-i)=-buse(8)
+        DD(1,2,-i)=-buse(8)
 !        DD(1,1,i)=0.0d0
 !        DD(2,2,i)=0.0d0
 !        DD(2,1,i)=0.0d0
 !        DD(1,2,i)=0.0d0
-        Dtilde(1,1,i)=b(6)
-        Dtilde(1,2,i)=b(8)
-        Dtilde(2,1,i)=-b(8)
-        Dtilde(2,2,i)=b(6)
-        Dtilde(1,1,-i)=b(6)
-        Dtilde(1,2,-i)=-b(8)
-        Dtilde(2,1,-i)=b(8)
-        Dtilde(2,2,-i)=b(6)
+        Dtilde(1,1,i)=buse(6)
+        Dtilde(1,2,i)=buse(8)
+        Dtilde(2,1,i)=-buse(8)
+        Dtilde(2,2,i)=buse(6)
+        Dtilde(1,1,-i)=buse(6)
+        Dtilde(1,2,-i)=-buse(8)
+        Dtilde(2,1,-i)=buse(8)
+        Dtilde(2,2,-i)=buse(6)
 
 !        Dtilde(1,1,i)=0.0d0
 !        Dtilde(1,2,i)=0.0d0
 !        Dtilde(2,1,i)=0.0d0
 !        Dtilde(2,2,i)=0.0d0
-        EE(1,1,i)= b(10)+b(11)
-        EE(2,2,i)=-b(10)+b(11)
-        EE(2,1,i)= b(12)-b(13)
-        EE(1,2,i)= b(12)+b(13)
-        EE(1,1,-i)= b(10)+b(11)
-        EE(2,2,-i)=-b(10)+b(11)
-        EE(2,1,-i)=-b(12)+b(13)
-        EE(1,2,-i)=-b(12)-b(13)
+        EE(1,1,i)= buse(10)+buse(11)
+        EE(2,2,i)=-buse(10)+buse(11)
+        EE(2,1,i)= buse(12)-buse(13)
+        EE(1,2,i)= buse(12)+buse(13)
+        EE(1,1,-i)= buse(10)+buse(11)
+        EE(2,2,-i)=-buse(10)+buse(11)
+        EE(2,1,-i)=-buse(12)+buse(13)
+        EE(1,2,-i)=-buse(12)-buse(13)
 
 !        ee(1,1,i)=1.0d0
 !        ee(2,2,i)=1.0d0
        sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
       write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
       enddo
+!-----------------READING SC BASE POTENTIALS-----------------------------
+      allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))      
+      allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
+      allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
+      allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
+      allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
+      allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
+      allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
+      allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
+      allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
+      allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
+      allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
+      allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
+      allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
+      allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
+      allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
+      allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
+
+
+      do i=1,ntyp_molec(1)
+       do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
+        write (*,*) "Im in ", i, " ", j
+        read(isidep_scbase,*) &
+        eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
+        chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
+         write(*,*) "eps",eps_scbase(i,j)
+        read(isidep_scbase,*) &
+       (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
+       chis_scbase(i,j,1),chis_scbase(i,j,2)
+        read(isidep_scbase,*) &
+       dhead_scbasei(i,j), &
+       dhead_scbasej(i,j), &
+       rborn_scbasei(i,j),rborn_scbasej(i,j)
+        read(isidep_scbase,*) &
+       (wdipdip_scbase(k,i,j),k=1,3), &
+       (wqdip_scbase(k,i,j),k=1,2)
+        read(isidep_scbase,*) &
+       alphapol_scbase(i,j), &
+       epsintab_scbase(i,j) 
+       END DO
+      END DO
+      allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
+      allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
+
+      do i=1,ntyp_molec(1)
+       do j=1,ntyp_molec(2)-1 
+          epsij=eps_scbase(i,j)
+          rrij=sigma_scbase(i,j)
+!          r0(i,j)=rrij
+!          r0(j,i)=rrij
+          rrij=rrij**expon
+!          epsij=eps(i,j)
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_scbase(i,j)=epsij*rrij*rrij
+          bb_scbase(i,j)=-sigeps*epsij*rrij
+        enddo
+       enddo
+!-----------------READING PEP BASE POTENTIALS-------------------
+      allocate(eps_pepbase(ntyp_molec(2)))
+      allocate(sigma_pepbase(ntyp_molec(2)))
+      allocate(chi_pepbase(ntyp_molec(2),2))
+      allocate(chipp_pepbase(ntyp_molec(2),2))
+      allocate(alphasur_pepbase(4,ntyp_molec(2)))
+      allocate(sigmap1_pepbase(ntyp_molec(2)))
+      allocate(sigmap2_pepbase(ntyp_molec(2)))
+      allocate(chis_pepbase(ntyp_molec(2),2))
+      allocate(wdipdip_pepbase(3,ntyp_molec(2)))
+
+
+       do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
+        write (*,*) "Im in ", i, " ", j
+        read(isidep_pepbase,*) &
+        eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
+        chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
+         write(*,*) "eps",eps_pepbase(j)
+        read(isidep_pepbase,*) &
+       (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
+       chis_pepbase(j,1),chis_pepbase(j,2)
+        read(isidep_pepbase,*) &
+       (wdipdip_pepbase(k,j),k=1,3)
+       END DO
+      allocate(aa_pepbase(ntyp_molec(2)))
+      allocate(bb_pepbase(ntyp_molec(2)))
+
+       do j=1,ntyp_molec(2)-1
+          epsij=eps_pepbase(j)
+          rrij=sigma_pepbase(j)
+!          r0(i,j)=rrij
+!          r0(j,i)=rrij
+          rrij=rrij**expon
+!          epsij=eps(i,j)
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_pepbase(j)=epsij*rrij*rrij
+          bb_pepbase(j)=-sigeps*epsij*rrij
+        enddo
+!--------------READING SC PHOSPHATE------------------------------------- 
+      allocate(eps_scpho(ntyp_molec(1)))
+      allocate(sigma_scpho(ntyp_molec(1)))
+      allocate(chi_scpho(ntyp_molec(1),2))
+      allocate(chipp_scpho(ntyp_molec(1),2))
+      allocate(alphasur_scpho(4,ntyp_molec(1)))
+      allocate(sigmap1_scpho(ntyp_molec(1)))
+      allocate(sigmap2_scpho(ntyp_molec(1)))
+      allocate(chis_scpho(ntyp_molec(1),2))
+      allocate(wqq_scpho(ntyp_molec(1)))
+      allocate(wqdip_scpho(2,ntyp_molec(1)))
+      allocate(alphapol_scpho(ntyp_molec(1)))
+      allocate(epsintab_scpho(ntyp_molec(1)))
+      allocate(dhead_scphoi(ntyp_molec(1)))
+      allocate(rborn_scphoi(ntyp_molec(1)))
+      allocate(rborn_scphoj(ntyp_molec(1)))
+
+!      j=1
+       do j=1,ntyp_molec(1) ! without U then we will take T for U
+        write (*,*) "Im in scpho ", i, " ", j
+        read(isidep_scpho,*) &
+        eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
+        chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
+         write(*,*) "eps",eps_scpho(j)
+        read(isidep_scpho,*) &
+       (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
+       chis_scpho(j,1),chis_scpho(j,2)
+        read(isidep_scpho,*) &
+       (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
+        read(isidep_scpho,*) &
+         epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j)
+       
+       END DO
+      allocate(aa_scpho(ntyp_molec(1)))
+      allocate(bb_scpho(ntyp_molec(1)))
+
+       do j=1,ntyp_molec(1)
+          epsij=eps_scpho(j)
+          rrij=sigma_scpho(j)
+!          r0(i,j)=rrij
+!          r0(j,i)=rrij
+          rrij=rrij**expon
+!          epsij=eps(i,j)
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_scpho(j)=epsij*rrij*rrij
+          bb_scpho(j)=-sigeps*epsij*rrij
+        enddo
+
+
+        read(isidep_peppho,*) &
+        eps_peppho,sigma_peppho
+        read(isidep_peppho,*) &
+       (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
+        read(isidep_peppho,*) &
+       (wqdip_peppho(k),k=1,2)
+
+          epsij=eps_peppho
+          rrij=sigma_peppho
+!          r0(i,j)=rrij
+!          r0(j,i)=rrij
+          rrij=rrij**expon
+!          epsij=eps(i,j)
+          sigeps=dsign(1.0D0,epsij)
+          epsij=dabs(epsij)
+          aa_peppho=epsij*rrij*rrij
+          bb_peppho=-sigeps*epsij*rrij
+
 
       allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
       bad(:,:)=0.0D0
 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
-           if (unres_pdb) then
+!           if (unres_pdb) then
 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
 !            print *,i,'tu dochodze'
-            call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
-            if (fail) then
-              e2(1)=0.0d0
-              e2(2)=1.0d0
-              e2(3)=0.0d0
-            endif !fail
+!            call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
+!            if (fail) then
+!              e2(1)=0.0d0
+!              e2(2)=1.0d0
+!              e2(3)=0.0d0
+!            endif !fail
 !            print *,i,'a tu?'
-            do j=1,3
-             c(j,i)=c(j,i-1)-1.9d0*e2(j)
-            enddo
-           else   !unres_pdb
+!            do j=1,3
+!             c(j,i)=c(j,i-1)-1.9d0*e2(j)
+!            enddo
+!           else   !unres_pdb
            do j=1,3
              dcj=(c(j,i-2)-c(j,i-3))/2.0
             if (dcj.eq.0) dcj=1.23591524223
              c(j,i)=c(j,i-1)+dcj
              c(j,nres+i)=c(j,i)
            enddo
-          endif   !unres_pdb
+!          endif   !unres_pdb
          else     !itype(i+1,1).eq.ntyp1
-          if (unres_pdb) then
+!          if (unres_pdb) then
 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
-            call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
-            if (fail) then
-              e2(1)=0.0d0
-              e2(2)=1.0d0
-              e2(3)=0.0d0
-            endif
+!            call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
+!            if (fail) then
+!              e2(1)=0.0d0
+!              e2(2)=1.0d0
+!              e2(3)=0.0d0
+!            endif
             do j=1,3
               c(j,i)=c(j,i+1)-1.9d0*e2(j)
             enddo
-          else !unres_pdb
+!          else !unres_pdb
            do j=1,3
             dcj=(c(j,i+3)-c(j,i+2))/2.0
             if (dcj.eq.0) dcj=1.23591524223
             c(j,i)=c(j,i+1)-dcj
             c(j,nres+i)=c(j,i)
            enddo
-          endif !unres_pdb
+!          endif !unres_pdb
          endif !itype(i+1,1).eq.ntyp1
         endif  !itype.eq.ntyp1
 
           if (molecule.eq.5) molecule=molecprev
         itype(nres,molecule)=ntyp1_molec(molecule)
         nres_molec(molecule)=nres_molec(molecule)+1
-        if (unres_pdb) then
+!        if (unres_pdb) then
 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
-          call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
-          if (fail) then
-            e2(1)=0.0d0
-            e2(2)=1.0d0
-            e2(3)=0.0d0
-          endif
-          do j=1,3
-            c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
-          enddo
-        else
+!          call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
+!          if (fail) then
+!            e2(1)=0.0d0
+!            e2(2)=1.0d0
+!            e2(3)=0.0d0
+!          endif
+!          do j=1,3
+!            c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
+!          enddo
+!        else
         do j=1,3
           dcj=(c(j,nres-2)-c(j,nres-3))/2.0
           c(j,nres)=c(j,nres-1)+dcj
           c(j,2*nres)=c(j,nres)
         enddo
-        endif
+!        endif
       endif
 !     print *,'I have',nres, nres_molec(:)
 
 ! NOW LETS ROCK! SORTING
       allocate(c_temporary(3,2*nres))
       allocate(itype_temporary(nres,5))
-      allocate(molnum(nres))
+      allocate(molnum(nres+1))
       allocate(istype_temp(nres))
        itype_temporary(:,:)=0
       seqalingbegin=1
       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('SIDEPAR_SCBASE',sidename_scbase)
+      open (isidep_scbase,file=sidename_scbase,status='old',action='read')
+      call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
+      open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
+      call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
+      open (isidep_scpho,file=pepname_scpho,status='old',action='read')
+      call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
+      open (isidep_peppho,file=pepname_peppho,status='old',action='read')
+
 
       call getenv_loc('LIPTRANPAR',liptranname)
       open (iliptranpar,file=liptranname,status='old',action='read')