--- /dev/null
+ subroutine minimize(etot,x,iretcode,nfun)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
+*********************************************************************
+* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
+* the calling subprogram. *
+* when d(i)=1.0, then v(35) is the length of the initial step, *
+* calculated in the usual pythagorean way. *
+* absolute convergence occurs when the function is within v(31) of *
+* zero. unless you know the minimum value in advance, abs convg *
+* is probably not useful. *
+* relative convergence is when the model predicts that the function *
+* will decrease by less than v(32)*abs(fun). *
+*********************************************************************
+ include 'COMMON.IOUNITS'
+ include 'COMMON.VAR'
+ include 'COMMON.GEO'
+ include 'COMMON.MINIM'
+ common /srutu/ icall
+ dimension iv(liv)
+ double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
+ double precision energia(0:n_ene)
+ external func,gradient,fdum
+ external func_restr,grad_restr
+ logical not_done,change,reduce
+c common /przechowalnia/ v
+
+ icall = 1
+
+ NOT_DONE=.TRUE.
+
+c DO WHILE (NOT_DONE)
+
+ call deflt(2,iv,liv,lv,v)
+* 12 means fresh start, dont call deflt
+ iv(1)=12
+* max num of fun calls
+ if (maxfun.eq.0) maxfun=500
+ iv(17)=maxfun
+* max num of iterations
+ if (maxmin.eq.0) maxmin=1000
+ iv(18)=maxmin
+* controls output
+ iv(19)=2
+* selects output unit
+ iv(21)=0
+ if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
+* 1 means to print out result
+ iv(22)=print_min_res
+* 1 means to print out summary stats
+ iv(23)=print_min_stat
+* 1 means to print initial x and d
+ iv(24)=print_min_ini
+* min val for v(radfac) default is 0.1
+ v(24)=0.1D0
+* max val for v(radfac) default is 4.0
+ v(25)=2.0D0
+c v(25)=4.0D0
+* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
+* the sumsl default is 0.1
+ v(26)=0.1D0
+* false conv if (act fnctn decrease) .lt. v(34)
+* the sumsl default is 100*machep
+ v(34)=v(34)/100.0D0
+* absolute convergence
+ if (tolf.eq.0.0D0) tolf=1.0D-4
+ v(31)=tolf
+* relative convergence
+ if (rtolf.eq.0.0D0) rtolf=1.0D-4
+ v(32)=rtolf
+* controls initial step size
+ v(35)=1.0D-1
+* large vals of d correspond to small components of step
+ do i=1,nphi
+ d(i)=1.0D-1
+ enddo
+ do i=nphi+1,nvar
+ d(i)=1.0D-1
+ enddo
+cd print *,'Calling SUMSL'
+c call var_to_geom(nvar,x)
+c call chainbuild
+c call etotal(energia(0))
+c etot = energia(0)
+ IF (mask_r) THEN
+ call x2xx(x,xx,nvar_restr)
+ call sumsl(nvar_restr,d,xx,func_restr,grad_restr,
+ & iv,liv,lv,v,idum,rdum,fdum)
+ call xx2x(x,xx)
+ ELSE
+ call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
+ ENDIF
+ etot=v(10)
+ iretcode=iv(1)
+cd print *,'Exit SUMSL; return code:',iretcode,' energy:',etot
+cd write (iout,'(/a,i4/)') 'SUMSL return code:',iv(1)
+c call intout
+c change=reduce(x)
+ call var_to_geom(nvar,x)
+c if (change) then
+c write (iout,'(a)') 'Reduction worked, minimizing again...'
+c else
+c not_done=.false.
+c endif
+ call chainbuild
+c call etotal(energia(0))
+c etot=energia(0)
+c call enerprint(energia(0))
+ nfun=iv(6)
+
+c write (*,*) 'Processor',MyID,' leaves MINIMIZE.'
+
+c ENDDO ! NOT_DONE
+
+ return
+ end
+#ifdef MPI
+c----------------------------------------------------------------------------
+ subroutine ergastulum
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include "mpif.h"
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.DERIV'
+ include 'COMMON.VAR'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.FFIELD'
+ include 'COMMON.INTERACT'
+ include 'COMMON.MD'
+ include 'COMMON.TIME1'
+ double precision z(maxres6),d_a_tmp(maxres6)
+ double precision edum(0:n_ene),time_order(0:10)
+ double precision Gcopy(maxres2,maxres2)
+ common /przechowalnia/ Gcopy
+ integer icall /0/
+C Workers wait for variables and NF, and NFL from the boss
+ iorder=0
+ do while (iorder.ge.0)
+c write (*,*) 'Processor',fg_rank,' CG group',kolor,
+c & ' receives order from Master'
+ time00=MPI_Wtime()
+ call MPI_Bcast(iorder,1,MPI_INTEGER,king,FG_COMM,IERR)
+ time_Bcast=time_Bcast+MPI_Wtime()-time00
+ if (icall.gt.4 .and. iorder.ge.0)
+ & time_order(iorder)=time_order(iorder)+MPI_Wtime()-time00
+ icall=icall+1
+c write (*,*)
+c & 'Processor',fg_rank,' completed receive MPI_BCAST order',iorder
+ if (iorder.eq.0) then
+ call zerograd
+ call etotal(edum)
+c write (2,*) "After etotal"
+c write (2,*) "dimen",dimen," dimen3",dimen3
+c call flush(2)
+ else if (iorder.eq.2) then
+ call zerograd
+ call etotal_short(edum)
+c write (2,*) "After etotal_short"
+c write (2,*) "dimen",dimen," dimen3",dimen3
+c call flush(2)
+ else if (iorder.eq.3) then
+ call zerograd
+ call etotal_long(edum)
+c write (2,*) "After etotal_long"
+c write (2,*) "dimen",dimen," dimen3",dimen3
+c call flush(2)
+ else if (iorder.eq.1) then
+ call sum_gradient
+c write (2,*) "After sum_gradient"
+c write (2,*) "dimen",dimen," dimen3",dimen3
+c call flush(2)
+ else if (iorder.eq.4) then
+ call ginv_mult(z,d_a_tmp)
+ else if (iorder.eq.5) then
+c Setup MD things for a slave
+ dimen=(nct-nnt+1)+nside
+ dimen1=(nct-nnt)+(nct-nnt+1)
+ dimen3=dimen*3
+c write (2,*) "dimen",dimen," dimen3",dimen3
+c call flush(2)
+ call int_bounds(dimen,igmult_start,igmult_end)
+ igmult_start=igmult_start-1
+ call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,
+ & ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
+ my_ng_count=igmult_end-igmult_start
+ call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,
+ & MPI_INTEGER,FG_COMM,IERROR)
+c write (2,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
+c write (2,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
+ myginv_ng_count=maxres2*my_ng_count
+c write (2,*) "igmult_start",igmult_start," igmult_end",
+c & igmult_end," my_ng_count",my_ng_count
+c call flush(2)
+ call MPI_Allgather(maxres2*igmult_start,1,MPI_INTEGER,
+ & nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
+ call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,
+ & nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
+c write (2,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
+c write (2,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
+c call flush(2)
+c call MPI_Barrier(FG_COMM,IERROR)
+ time00=MPI_Wtime()
+ call MPI_Scatterv(ginv(1,1),nginv_counts(0),
+ & nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
+ & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
+#ifdef TIMING
+ time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
+#endif
+ do i=1,dimen
+ do j=1,2*my_ng_count
+ ginv(j,i)=gcopy(i,j)
+ enddo
+ enddo
+c write (2,*) "dimen",dimen," dimen3",dimen3
+c write (2,*) "End MD setup"
+c call flush(2)
+c write (iout,*) "My chunk of ginv_block"
+c call MATOUT2(my_ng_count,dimen3,maxres2,maxers2,ginv_block)
+ else if (iorder.eq.6) then
+ call int_from_cart1(.false.)
+ else if (iorder.eq.7) then
+ call chainbuild_cart
+ else if (iorder.eq.8) then
+ call intcartderiv
+ else if (iorder.eq.9) then
+ call fricmat_mult(z,d_a_tmp)
+ else if (iorder.eq.10) then
+ call setup_fricmat
+ endif
+ enddo
+ write (*,*) 'Processor',fg_rank,' CG group',kolor,
+ & ' absolute rank',myrank,' leves ERGASTULUM.'
+ write(*,*)'Processor',fg_rank,' wait times for respective orders',
+ & (' order[',i,']',time_order(i),i=0,10)
+ return
+ end
+#endif
+************************************************************************
+ subroutine func(n,x,nf,f,uiparm,urparm,ufparm)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.DERIV'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ common /chuju/ jjj
+ double precision energia(0:n_ene)
+ integer jjj
+ double precision ufparm
+ external ufparm
+ integer uiparm(1)
+ real*8 urparm(1)
+ dimension x(maxvar)
+c if (jjj.gt.0) then
+c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
+c endif
+ nfl=nf
+ icg=mod(nf,2)+1
+cd print *,'func',nf,nfl,icg
+ call var_to_geom(n,x)
+ call zerograd
+ call chainbuild
+cd write (iout,*) 'ETOTAL called from FUNC'
+ call etotal(energia(0))
+ call sum_gradient
+ f=energia(0)
+c if (jjj.gt.0) then
+c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
+c write (iout,*) 'f=',etot
+c jjj=0
+c endif
+ return
+ end
+************************************************************************
+ subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.DERIV'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ common /chuju/ jjj
+ double precision energia(0:n_ene)
+ integer jjj
+ double precision ufparm
+ external ufparm
+ integer uiparm(1)
+ real*8 urparm(1)
+ dimension x(maxvar)
+c if (jjj.gt.0) then
+c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
+c endif
+ nfl=nf
+ icg=mod(nf,2)+1
+ call var_to_geom_restr(n,x)
+ call zerograd
+ call chainbuild
+cd write (iout,*) 'ETOTAL called from FUNC'
+ call etotal(energia(0))
+ call sum_gradient
+ f=energia(0)
+c if (jjj.gt.0) then
+c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
+c write (iout,*) 'f=',etot
+c jjj=0
+c endif
+ return
+ end
+c-------------------------------------------------------
+ subroutine x2xx(x,xx,n)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.VAR'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ double precision xx(maxvar),x(maxvar)
+
+ do i=1,nvar
+ varall(i)=x(i)
+ enddo
+
+ ig=0
+ igall=0
+ do i=4,nres
+ igall=igall+1
+ if (mask_phi(i).eq.1) then
+ ig=ig+1
+ xx(ig)=x(igall)
+ endif
+ enddo
+
+ do i=3,nres
+ igall=igall+1
+ if (mask_theta(i).eq.1) then
+ ig=ig+1
+ xx(ig)=x(igall)
+ endif
+ enddo
+
+ do ij=1,2
+ do i=2,nres-1
+ if (itype(i).ne.10) then
+ igall=igall+1
+ if (mask_side(i).eq.1) then
+ ig=ig+1
+ xx(ig)=x(igall)
+ endif
+ endif
+ enddo
+ enddo
+
+ n=ig
+
+ return
+ end
+
+ subroutine xx2x(x,xx)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.VAR'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ double precision xx(maxvar),x(maxvar)
+
+ do i=1,nvar
+ x(i)=varall(i)
+ enddo
+
+ ig=0
+ igall=0
+ do i=4,nres
+ igall=igall+1
+ if (mask_phi(i).eq.1) then
+ ig=ig+1
+ x(igall)=xx(ig)
+ endif
+ enddo
+
+ do i=3,nres
+ igall=igall+1
+ if (mask_theta(i).eq.1) then
+ ig=ig+1
+ x(igall)=xx(ig)
+ endif
+ enddo
+
+ do ij=1,2
+ do i=2,nres-1
+ if (itype(i).ne.10) then
+ igall=igall+1
+ if (mask_side(i).eq.1) then
+ ig=ig+1
+ x(igall)=xx(ig)
+ endif
+ endif
+ enddo
+ enddo
+
+ return
+ end
+
+c----------------------------------------------------------
+ subroutine minim_dc(etot,iretcode,nfun)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.VAR'
+ include 'COMMON.GEO'
+ include 'COMMON.MINIM'
+ include 'COMMON.CHAIN'
+ dimension iv(liv)
+ double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
+c common /przechowalnia/ v
+
+ double precision energia(0:n_ene)
+ external func_dc,grad_dc,fdum
+ logical not_done,change,reduce
+ double precision g(maxvar),f1
+
+ call deflt(2,iv,liv,lv,v)
+* 12 means fresh start, dont call deflt
+ iv(1)=12
+* max num of fun calls
+ if (maxfun.eq.0) maxfun=500
+ iv(17)=maxfun
+* max num of iterations
+ if (maxmin.eq.0) maxmin=1000
+ iv(18)=maxmin
+* controls output
+ iv(19)=2
+* selects output unit
+ iv(21)=0
+ if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
+* 1 means to print out result
+ iv(22)=print_min_res
+* 1 means to print out summary stats
+ iv(23)=print_min_stat
+* 1 means to print initial x and d
+ iv(24)=print_min_ini
+* min val for v(radfac) default is 0.1
+ v(24)=0.1D0
+* max val for v(radfac) default is 4.0
+ v(25)=2.0D0
+c v(25)=4.0D0
+* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
+* the sumsl default is 0.1
+ v(26)=0.1D0
+* false conv if (act fnctn decrease) .lt. v(34)
+* the sumsl default is 100*machep
+ v(34)=v(34)/100.0D0
+* absolute convergence
+ if (tolf.eq.0.0D0) tolf=1.0D-4
+ v(31)=tolf
+* relative convergence
+ if (rtolf.eq.0.0D0) rtolf=1.0D-4
+ v(32)=rtolf
+* controls initial step size
+ v(35)=1.0D-1
+* large vals of d correspond to small components of step
+ do i=1,6*nres
+ d(i)=1.0D-1
+ enddo
+
+ k=0
+ do i=1,nres-1
+ do j=1,3
+ k=k+1
+ x(k)=dc(j,i)
+ enddo
+ enddo
+ do i=2,nres-1
+ if (ialph(i,1).gt.0) then
+ do j=1,3
+ k=k+1
+ x(k)=dc(j,i+nres)
+ enddo
+ endif
+ enddo
+
+ call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
+
+ k=0
+ do i=1,nres-1
+ do j=1,3
+ k=k+1
+ dc(j,i)=x(k)
+ enddo
+ enddo
+ do i=2,nres-1
+ if (ialph(i,1).gt.0) then
+ do j=1,3
+ k=k+1
+ dc(j,i+nres)=x(k)
+ enddo
+ endif
+ enddo
+ call chainbuild_cart
+
+cd call zerograd
+cd nf=0
+cd call func_dc(k,x,nf,f,idum,rdum,fdum)
+cd call grad_dc(k,x,nf,g,idum,rdum,fdum)
+cd
+cd do i=1,k
+cd x(i)=x(i)+1.0D-5
+cd call func_dc(k,x,nf,f1,idum,rdum,fdum)
+cd x(i)=x(i)-1.0D-5
+cd print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
+cd enddo
+
+ etot=v(10)
+ iretcode=iv(1)
+ nfun=iv(6)
+ return
+ end
+
+ subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.DERIV'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.CHAIN'
+ include 'COMMON.VAR'
+ double precision energia(0:n_ene)
+ double precision ufparm
+ external ufparm
+ integer uiparm(1)
+ real*8 urparm(1)
+ dimension x(maxvar)
+ nfl=nf
+cbad icg=mod(nf,2)+1
+ icg=1
+
+ k=0
+ do i=1,nres-1
+ do j=1,3
+ k=k+1
+ dc(j,i)=x(k)
+ enddo
+ enddo
+ do i=2,nres-1
+ if (ialph(i,1).gt.0) then
+ do j=1,3
+ k=k+1
+ dc(j,i+nres)=x(k)
+ enddo
+ endif
+ enddo
+ call chainbuild_cart
+
+ call zerograd
+ call etotal(energia(0))
+ f=energia(0)
+
+cd print *,'func_dc ',nf,nfl,f
+
+ return
+ end
+
+ subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ include 'COMMON.SETUP'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.VAR'
+ include 'COMMON.INTERACT'
+ include 'COMMON.FFIELD'
+ include 'COMMON.MD'
+ include 'COMMON.IOUNITS'
+ external ufparm
+ integer uiparm(1),k
+ double precision urparm(1)
+ dimension x(maxvar),g(maxvar)
+c
+c
+c
+cbad icg=mod(nf,2)+1
+ icg=1
+cd print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
+ if (nf-nfl+1) 20,30,40
+ 20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
+cd print *,20
+ if (nf.eq.0) return
+ goto 40
+ 30 continue
+cd print *,30
+ k=0
+ do i=1,nres-1
+ do j=1,3
+ k=k+1
+ dc(j,i)=x(k)
+ enddo
+ enddo
+ do i=2,nres-1
+ if (ialph(i,1).gt.0) then
+ do j=1,3
+ k=k+1
+ dc(j,i+nres)=x(k)
+ enddo
+ endif
+ enddo
+ call chainbuild_cart
+
+C
+C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
+C
+ 40 call cartgrad
+cd print *,40
+ k=0
+ do i=1,nres-1
+ do j=1,3
+ k=k+1
+ g(k)=gcart(j,i)
+ enddo
+ enddo
+ do i=2,nres-1
+ if (ialph(i,1).gt.0) then
+ do j=1,3
+ k=k+1
+ g(k)=gxcart(j,i)
+ enddo
+ endif
+ enddo
+
+ return
+ end
projects / unres.git / blobdiff