corrections

[unres.git] / source / unres / src-HCD-5D / contact_cp2.F

diff --git a/source/unres/src-HCD-5D/contact_cp2.F b/source/unres/src-HCD-5D/contact_cp2.F
new file mode 100644 (file)
index 0000000..785c8cb
--- /dev/null
+++ b/source/unres/src-HCD-5D/contact_cp2.F
@@ -0,0 +1,148 @@
+      subroutine contact_cp2(var,var2,iff,ieval,in_pdb)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.FFIELD'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.DISTFIT'
+      include 'COMMON.VAR'
+      include 'COMMON.CHAIN'
+      include 'COMMON.MINIM'
+
+      character*50 linia
+      integer nf,ij(4)
+      double precision var(maxvar),var2(maxvar)
+      double precision time0,time1
+      integer iff(maxres),ieval      
+      double precision theta1(maxres),phi1(maxres),alph1(maxres),     
+     &                 omeg1(maxres)                             
+      
+
+       call var_to_geom(nvar,var)
+       call chainbuild                                                           
+       nhpb0=nhpb
+       ind=0                                                                     
+       do i=1,nres-3                                                             
+         do j=i+3,nres                                                           
+          ind=ind+1                                                              
+          if ( iff(i).eq.1.and.iff(j).eq.1 ) then                                           
+            d0(ind)=DIST(i,j)                                                     
+            w(ind)=10.0                                                           
+            nhpb=nhpb+1                                                           
+            ihpb(nhpb)=i                                                          
+            jhpb(nhpb)=j                                                          
+            forcon(nhpb)=10.0                                                     
+            dhpb(nhpb)=d0(ind)                                                    
+          else
+            w(ind)=0.0
+          endif                                                                  
+         enddo                                                                   
+       enddo                                    
+       call hpb_partition
+
+       do i=1,nres                                                               
+        theta1(i)=theta(i)                                                      
+        phi1(i)=phi(i)                                                          
+        alph1(i)=alph(i)                                                        
+        omeg1(i)=omeg(i)                                                        
+       enddo                      
+
+       call var_to_geom(nvar,var2)
+
+       do i=1,nres                                                             
+          if ( iff(i).eq.1 ) then                                           
+              theta(i)=theta1(i)                                                      
+              phi(i)=phi1(i)                                                          
+              alph(i)=alph1(i)                                                        
+              omeg(i)=omeg1(i)                       
+          endif
+       enddo
+
+       call chainbuild
+cd       call write_pdb(3,'combined structure',0d0)
+cd       time0=MPI_WTIME()
+       
+       NX=NRES-3                                                                 
+       NY=((NRES-4)*(NRES-5))/2 
+       call distfit(.true.,200)
+   
+cd       time1=MPI_WTIME()
+cd       write (iout,'(a,f6.2,a)') ' Time for distfit ',time1-time0,' sec'
+
+       ipot0=ipot
+       maxmin0=maxmin
+       maxfun0=maxfun
+       wstrain0=wstrain
+
+       ipot=6
+       maxmin=2000
+       maxfun=5000
+       call geom_to_var(nvar,var)
+cd       time0=MPI_WTIME()
+       call minimize(etot,var,iretcode,nfun)                               
+       write(iout,*)'SUMSL return code is',iretcode,' eval SOFT',nfun   
+
+cd       time1=MPI_WTIME()
+cd       write (iout,'(a,f6.2,f8.2,a)')' Time for soft min.',time1-time0,
+cd     &         nfun/(time1-time0),' SOFT eval/s'
+        call var_to_geom(nvar,var)
+        call chainbuild
+
+
+        iwsk=0
+        nf=0
+        if (iff(1).eq.1) then
+          iwsk=1
+          nf=nf+1
+          ij(nf)=0
+        endif
+        do i=2,nres
+           if ( iwsk.eq.0.and.iff(i-1).eq.0.and.iff(i).eq.1 ) then
+             iwsk=1
+             nf=nf+1
+             ij(nf)=i
+           endif
+           if ( iwsk.eq.1.and.iff(i-1).eq.1.and.iff(i).eq.0 ) then
+             iwsk=0
+             nf=nf+1
+             ij(nf)=i-1
+           endif
+        enddo
+        if (iff(nres).eq.1) then
+          nf=nf+1
+          ij(nf)=nres
+        endif
+
+
+cd        write(linia,'(a6,i3,a1,i3,a1,i3,a1,i3)')
+cd     &                     "select",ij(1),"-",ij(2),
+cd     &                         ",",ij(3),"-",ij(4)
+cd        call write_pdb(in_pdb,linia,etot)
+
+
+       ipot=ipot0
+       maxmin=maxmin0
+       maxfun=maxfun0
+cd       time0=MPI_WTIME()
+       call minimize(etot,var,iretcode,nfun)
+cd       write(iout,*)'SUMSL DIST return code is',iretcode,' eval ',nfun
+       ieval=nfun
+
+cd       time1=MPI_WTIME()
+cd       write (iout,'(a,f6.2,f8.2,a)')' Time for DIST min.',time1-time0,
+cd     &         nfun/(time1-time0),' eval/s'
+cd       call var_to_geom(nvar,var)
+cd       call chainbuild
+cd       call write_pdb(6,'dist structure',etot)
+
+
+       nhpb= nhpb0                                                                  
+       link_start=1                                                            
+       link_end=nhpb     
+       wstrain=wstrain0
+
+       return
+       end