ENERGY_DEC printout works for ebend in E0LL2Y forcefield
[unres.git] / source / unres / src_MD / readrtns.F
index 4fd5a31..887e89f 100644 (file)
@@ -2719,8 +2719,10 @@ c    &    sigma_odl_temp(maxres,maxres,max_template)
       character*2 kic2
       character*24 model_ki_dist, model_ki_angle
       character*500 controlcard
-      integer ki, i, j, k, l, ii_in_use(maxdim)
+      integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
       logical lprn /.true./
+      integer ilen
+      external ilen
 c
 c     FP - Nov. 2014 Temporary specifications for new vars
 c
@@ -2772,6 +2774,14 @@ c
       lim_theta=0
       lim_xx=0
 c
+      write(iout,*) 'nnt=',nnt,'nct=',nct
+c
+      do i = nnt,nct
+        do k=1,constr_homology
+         idomain(k,i)=0
+        enddo
+      enddo
+
       ii=0
       do i = nnt,nct-2 
         do j=i+2,nct 
 
       do k=1,constr_homology
 
-c
-c Files containing res sim or local scores (former containing sigmas)
-c
-
-        write(kic2,'(bz,i2.2)') k
-
-        tpl_k_rescore="template"//kic2//".sco"
-
-c
-c From read_dist_constr (commented out 25/11/2014 <-> res sim)
-c
-c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
-          open (ientin,file=tpl_k_rescore,status='old')
-          do irec=1,maxdim ! loop for reading res sim 
-            if (irec.eq.1) then
-               rescore(k,irec)=0.0d0
-               goto 1301 
-            endif
-            if (read2sigma) then
-             read (ientin,*,end=1401) rescore2_tmp,rescore_tmp,
-     &                                idomain(k,irec)
-             rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0)
-             rescore2(k,irec)=0.5d0*(rescore2_tmp+0.5d0)
-            else
-             idomain(k,irec)=1
-             read (ientin,*,end=1401) rescore_tmp
-
-c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
-             rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
-c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
-            endif
- 1301     continue
-          enddo  
- 1401   continue
-          close (ientin)        
-
         read(inp,'(a)') pdbfile
 c  Next stament causes error upon compilation (?)
 c       if(me.eq.king.or. .not. out1file)
 c         write (iout,'(2a)') 'PDB data will be read from file ',
 c    &   pdbfile(:ilen(pdbfile))
+         write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
+     &  pdbfile(:ilen(pdbfile))
         open(ipdbin,file=pdbfile,status='old',err=33)
         goto 34
-  33    write (iout,'(a)') 'Error opening PDB file.'
+  33    write (iout,'(a,5x,a)') 'Error opening PDB file',
+     &  pdbfile(:ilen(pdbfile))
         stop
   34    continue
 c        print *,'Begin reading pdb data'
+c
+c Files containing res sim or local scores (former containing sigmas)
+c
+
+        write(kic2,'(bz,i2.2)') k
+
+        tpl_k_rescore="template"//kic2//".sco"
 
         unres_pdb=.false.
         if (read2sigma) then
-         call readpdb_template(k)
+          call readpdb_template(k)
         else
-         call readpdb
+          call readpdb
         endif
 c
+c     Distance restraints
+c
+c          ... --> odl(k,ii)
 C Copy the coordinates from reference coordinates (?)
         do i=1,2*nres
           do j=1,3
@@ -2844,17 +2831,35 @@ C Copy the coordinates from reference coordinates (?)
 c           write (iout,*) "c(",j,i,") =",c(j,i)
           enddo
         enddo
+c
+c From read_dist_constr (commented out 25/11/2014 <-> res sim)
+c
+c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
+          open (ientin,file=tpl_k_rescore,status='old')
+          if (nnt.gt.1) rescore(k,1)=0.0d0
+          do irec=nnt,maxdim ! loop for reading res sim 
+            if (read2sigma) then
+             read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
+     &                                idomain_tmp
+             i_tmp=i_tmp+nnt-1
+             idomain(k,i_tmp)=idomain_tmp
+             rescore(k,i_tmp)=0.5d0*(rescore_tmp+0.5d0)
+             rescore2(k,i_tmp)=0.5d0*(rescore2_tmp+0.5d0)
+            else
+             idomain(k,irec)=1
+             read (ientin,*,end=1401) rescore_tmp
 
-
-
+c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
+             rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
+c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
+            endif
+          enddo  
+ 1401   continue
+        close (ientin)        
         if (waga_dist.ne.0.0d0) then
           ii=0
-          do i = nnt,nct-2 ! right? without parallel.
-            do j=i+2,nct ! right?
-c         do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology 
-c           do j=i+2,nres ! ibid
-c         do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
-c           do j=i+2,nct ! ibid
+          do i = nnt,nct-2 
+            do j=i+2,nct 
 
             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0) then
                   
               z12=c(3,i)-c(3,j)
               distal=dsqrt(x12*x12+y12*y12+z12*z12)
               odl(k,ii)=distal
-c
-c             odl(k,ii)=dist(i,j)
-c             write (iout,*) "dist(",i,j,") =",dist(i,j)
-c             write (iout,*) "distal = ",distal
-c             write (iout,*) "odl(",k,ii,") =",odl(k,ii)
-c            write(iout,*) "rescore(",k,i,") =",rescore(k,i),
-c    &                      "rescore(",k,j,") =",rescore(k,j)
-c
-c  Calculation of sigma from res sim
-c
-c             if (odl(k,ii).le.6.0d0) then
-c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
-c    Other functional forms possible depending on odl(k,ii), eg.
-c
-           if (read2sigma) then
-            sigma_odl(k,ii)=0
-            do ik=i,j
-              sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
-            enddo
-            sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
-            if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = 
-     &         sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
-           else
-            if (odl(k,ii).le.dist_cut) then              
-              sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
-c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
-            else
+              if (read2sigma) then
+                sigma_odl(k,ii)=0
+                do ik=i,j
+                 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
+                enddo
+                sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
+                if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = 
+     &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
+              else
+                if (odl(k,ii).le.dist_cut) then
+                 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
+c                sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
+                else
 #ifdef OLDSIGMA
-              sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error 
+                 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error 
      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
 #else
-              sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error 
+                 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error 
      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
 #endif
-            endif
-           endif
+                endif
+              endif
 c   Following expr replaced by a positive exp argument
 c             sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
 
 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
-            
 c
               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
 c             sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
@@ -2945,10 +2935,7 @@ c    &                            sigma_dih(k,i+nnt-1)
 c         enddo
 c1402   continue
 c         close (ientin)
-          do i = nnt+3,nct ! right? without parallel.
-c         do i=1,nres ! alternative for bounds acc to readpdb?
-c         do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
-c         do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
+          do i = nnt+3,nct 
             if (idomain(k,i).eq.0) then 
                sigma_dih(k,i)=0.0
                cycle
@@ -2962,8 +2949,8 @@ c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
 
             sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
-     &                     rescore(k,i-2)+rescore(k,i-3)  !  right expression ?
-c
+     &                     rescore(k,i-2)+rescore(k,i-3)
+            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
@@ -2999,7 +2986,8 @@ c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
              sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
-     &                        rescore(k,i-2) !  right expression ?
+     &                        rescore(k,i-2)
+             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
 
 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
@@ -3017,7 +3005,6 @@ c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use
 c    &                          sigma_d(k,i+nnt-1)
 c         enddo
 c1404   continue
-          
 
           do i = nnt,nct ! right? without parallel.
 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
@@ -3040,6 +3027,7 @@ c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
 c              read (ientin,*) sigma_d(k,i) ! 1st variant
+               if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
           enddo
           lim_xx=nct-nnt+1 
         endif