critical bugfix for nucleic sequence reading

[unres4.git] / source / unres / io.F90

diff --git a/source/unres/io.F90 b/source/unres/io.F90
index cc8b681..b2ebc22 100644 (file)
--- a/source/unres/io.F90
+++ b/source/unres/io.F90
@@ -736,10 +736,12 @@
 
       real(kind=8),dimension(3,maxres2+2) :: c_alloc
       real(kind=8),dimension(3,0:maxres2) :: dc_alloc
 
       real(kind=8),dimension(3,maxres2+2) :: c_alloc
       real(kind=8),dimension(3,0:maxres2) :: dc_alloc

+      real(kind=8),dimension(:,:), allocatable :: secprob

       integer,dimension(maxres) :: itype_alloc
 
       integer :: iti,nsi,maxsi,itrial,itmp
       integer,dimension(maxres) :: itype_alloc
 
       integer :: iti,nsi,maxsi,itrial,itmp

-      real(kind=8) :: wlong,scalscp,co,ssscale
+      real(kind=8) :: wlong,scalscp,co,ssscale,phihel,phibet,sigmahel,&
+      sigmabet,sumv

       allocate(weights(n_ene))
 !-----------------------------
       allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres
       allocate(weights(n_ene))
 !-----------------------------
       allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres

@@ -782,9 +784,10 @@
       call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0)
       call reada(weightcard,'WSBLOC',wsbloc,0.0D0)
       call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0)
       call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0)
       call reada(weightcard,'WSBLOC',wsbloc,0.0D0)
       call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0)

+!      print *,"KUR..",wtor_nucl

       call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0)
       call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0)
       call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0)
       call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0)

-      call reada(weightcard,'WCORR3_NUCL',wcorr3_nucl,0.0D0)
+      call reada(weightcard,'WCORR3_NUC',wcorr3_nucl,0.0D0)

       call reada(weightcard,'WBOND',wbond,1.0D0)
       call reada(weightcard,'WTOR',wtor,1.0D0)
       call reada(weightcard,'WTORD',wtor_d,1.0D0)
       call reada(weightcard,'WBOND',wbond,1.0D0)
       call reada(weightcard,'WTOR',wtor,1.0D0)
       call reada(weightcard,'WTORD',wtor_d,1.0D0)

@@ -830,6 +833,26 @@
       weights(17)=wbond
       weights(18)=scal14
       weights(21)=wsccor
       weights(17)=wbond
       weights(18)=scal14
       weights(21)=wsccor

+          weights(26)=wvdwpp_nucl
+          weights(27)=welpp
+          weights(28)=wvdwpsb
+          weights(29)=welpsb
+          weights(30)=wvdwsb
+          weights(31)=welsb
+          weights(32)=wbond_nucl
+          weights(33)=wang_nucl
+          weights(34)=wsbloc
+          weights(35)=wtor_nucl
+          weights(36)=wtor_d_nucl
+          weights(37)=wcorr_nucl
+          weights(38)=wcorr3_nucl
+          weights(41)=wcatcat
+          weights(42)=wcatprot
+          weights(46)=wscbase
+          weights(47)=wpepbase
+          weights(48)=wscpho
+          weights(49)=wpeppho
+

       if(me.eq.king.or..not.out1file) &
        write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
       if(me.eq.king.or..not.out1file) &
        write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
         wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&

@@ -953,6 +976,8 @@
   34    continue
 !        print *,'Begin reading pdb data'
         call readpdb
   34    continue
 !        print *,'Begin reading pdb data'
         call readpdb

+!        call int_from_cart1(.true.)
+

 !        print *,'Finished reading pdb data'
         if(me.eq.king.or..not.out1file) &
          write (iout,'(a,i3,a,i3)')'nsup=',nsup,&
 !        print *,'Finished reading pdb data'
         if(me.eq.king.or..not.out1file) &
          write (iout,'(a,i3,a,i3)')'nsup=',nsup,&

@@ -1031,7 +1056,7 @@
         do i=1,nres_molec(molec)
           itmp=itmp+1
           istype(itmp)=sugarcode(sequence(i,molec)(1:1),i)
         do i=1,nres_molec(molec)
           itmp=itmp+1
           istype(itmp)=sugarcode(sequence(i,molec)(1:1),i)

-          itype(itmp,molec)=rescode(i,sequence(i,molec)(2:2),iscode,molec)
+          itype(itmp,molec)=rescode(i,sequence(i,molec)(1:2),iscode,molec)

         enddo
        endif
 
         enddo
        endif
 

@@ -1052,7 +1077,7 @@
        nres=0
        do i=1,5
         nres=nres+nres_molec(i)
        nres=0
        do i=1,5
         nres=nres+nres_molec(i)

-        print *,nres_molec(i)
+        print *,"nres_molec",nres,nres_molec(i)

        enddo
        
 ! Assign initial virtual bond lengths
        enddo
        
 ! Assign initial virtual bond lengths

@@ -1067,14 +1092,18 @@
          enddo
 !        print *,nres_molec(i)
         endif
          enddo
 !        print *,nres_molec(i)
         endif

-        if(.not.allocated(vbld)) allocate(vbld(2*nres))
+        print *,nres,"nres"
+        if(.not.allocated(vbld)) then
+           print *, "I DO ENTER" 
+           allocate(vbld(2*nres))
+        endif

         if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres))
         do i=2,nres
           vbld(i)=vbl
           vbld_inv(i)=vblinv
         enddo
         do i=2,nres-1
         if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres))
         do i=2,nres
           vbld(i)=vbl
           vbld_inv(i)=vblinv
         enddo
         do i=2,nres-1

-          print *, "molnum",molnum(i),itype(i,molnum(i)) 
+!          print *, "molnum",molnum(i),itype(i,molnum(i)),nres,i 

           vbld(i+nres)=dsc(iabs(itype(i,molnum(i))))
           vbld_inv(i+nres)=dsc_inv(iabs(itype(i,molnum(i))))
 !          write (iout,*) "i",i," itype",itype(i,1),
           vbld(i+nres)=dsc(iabs(itype(i,molnum(i))))
           vbld_inv(i+nres)=dsc_inv(iabs(itype(i,molnum(i))))
 !          write (iout,*) "i",i," itype",itype(i,1),

@@ -1143,7 +1172,7 @@
       endif
       call read_bridge
 !--------------------------------
       endif
       call read_bridge
 !--------------------------------

-       print *,"tu dochodze"
+!       print *,"tu dochodze"

 ! znamy nres oraz nss można zaalokowac potrzebne tablice
       call alloc_geo_arrays
       call alloc_ener_arrays
 ! znamy nres oraz nss można zaalokowac potrzebne tablice
       call alloc_geo_arrays
       call alloc_ener_arrays

@@ -1155,6 +1184,7 @@
       enddo
       read (inp,*) ndih_constr
       if (ndih_constr.gt.0) then
       enddo
       read (inp,*) ndih_constr
       if (ndih_constr.gt.0) then

+        raw_psipred=.false.

         allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr)
         allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr)
         allocate(ftors(ndih_constr)) !(maxdih_constr)
         allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr)
         allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr)
         allocate(ftors(ndih_constr)) !(maxdih_constr)

@@ -1181,6 +1211,45 @@
           phibound(1,ii) = phi0(i)-drange(i)
           phibound(2,ii) = phi0(i)+drange(i)
         enddo 
           phibound(1,ii) = phi0(i)-drange(i)
           phibound(2,ii) = phi0(i)+drange(i)
         enddo 

+      else if (ndih_constr.lt.0) then
+        raw_psipred=.true.
+        allocate(secprob(3,nres))
+        allocate(vpsipred(3,nres))
+        allocate(sdihed(2,nres))
+        call card_concat(weightcard,.true.)
+        call reada(weightcard,"PHIHEL",phihel,50.0D0)
+        call reada(weightcard,"PHIBET",phibet,180.0D0)
+        call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
+        call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
+        call reada(weightcard,"WDIHC",wdihc,0.591D0)
+        write (iout,*) "Weight of dihedral angle restraints",wdihc
+        read(inp,'(9x,3f7.3)') &
+          (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
+        write (iout,*) "The secprob array"
+        do i=nnt,nct
+          write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
+        enddo
+        ndih_constr=0
+        do i=nnt+3,nct
+          if (itype(i-3,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1 &
+          .and. itype(i-1,1).ne.ntyp1 .and. itype(i,1).ne.ntyp1) then
+            ndih_constr=ndih_constr+1
+            idih_constr(ndih_constr)=i
+            sumv=0.0d0
+            do j=1,3
+              vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
+              sumv=sumv+vpsipred(j,ndih_constr)
+            enddo
+            do j=1,3
+              vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
+            enddo
+            phibound(1,ndih_constr)=phihel*deg2rad
+            phibound(2,ndih_constr)=phibet*deg2rad
+            sdihed(1,ndih_constr)=sigmahel*deg2rad
+            sdihed(2,ndih_constr)=sigmabet*deg2rad
+          endif
+        enddo
+

       endif
       if (with_theta_constr) then
 !C with_theta_constr is keyword allowing for occurance of theta constrains
       endif
       if (with_theta_constr) then
 !C with_theta_constr is keyword allowing for occurance of theta constrains

@@ -1371,7 +1440,10 @@
             enddo
             return
           else
             enddo
             return
           else

-            call read_angles(inp,*36)
+           write(iout,*) "read angles from input" 
+           call read_angles(inp,*36)
+            call chainbuild
+

           endif
           goto 37
    36     write (iout,'(a)') 'Error reading angle file.'
           endif
           goto 37
    36     write (iout,'(a)') 'Error reading angle file.'

@@ -1508,6 +1580,10 @@
           enddo
       endif
       if (i2ndstr.gt.0) call secstrp2dihc
           enddo
       endif
       if (i2ndstr.gt.0) call secstrp2dihc

+      if (indpdb.gt.0) then 
+          write(iout,*) "WCHODZE TU!!"
+          call int_from_cart1(.true.)
+      endif

 !      call geom_to_var(nvar,x)
 !      call etotal(energia(0))
 !      call enerprint(energia(0))
 !      call geom_to_var(nvar,x)
 !      call etotal(energia(0))
 !      call enerprint(energia(0))