NEW CORR added to UNRES4

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

diff --git a/source/unres/io.F90 b/source/unres/io.F90
index cc8b681..abce524 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(:,:), allocatable :: secprob
       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
@@ -1155,6 +1157,7 @@
       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)
@@ -1181,6 +1184,45 @@
           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
@@ -1371,7 +1413,10 @@
             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.'