1 subroutine minimize(etot,x,iretcode,nfun)
2 implicit real*8 (a-h,o-z)
4 parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
5 *********************************************************************
6 * OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
7 * the calling subprogram. *
8 * when d(i)=1.0, then v(35) is the length of the initial step, *
9 * calculated in the usual pythagorean way. *
10 * absolute convergence occurs when the function is within v(31) of *
11 * zero. unless you know the minimum value in advance, abs convg *
12 * is probably not useful. *
13 * relative convergence is when the model predicts that the function *
14 * will decrease by less than v(32)*abs(fun). *
15 *********************************************************************
16 include 'COMMON.IOUNITS'
19 include 'COMMON.MINIM'
22 double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
23 double precision energia(0:n_ene)
24 external func,gradient,fdum
25 external func_restr,grad_restr
26 logical not_done,change,reduce
27 c common /przechowalnia/ v
30 c print *,"liv",liv," lv",lv
31 call deflt(2,iv,liv,lv,v)
32 * 12 means fresh start, dont call deflt
34 * max num of fun calls
35 if (maxfun.eq.0) maxfun=500
37 * max num of iterations
38 if (maxmin.eq.0) maxmin=1000
44 if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
45 * 1 means to print out result
47 * 1 means to print out summary stats
49 * 1 means to print initial x and d
51 * min val for v(radfac) default is 0.1
53 * max val for v(radfac) default is 4.0
56 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
57 * the sumsl default is 0.1
59 * false conv if (act fnctn decrease) .lt. v(34)
60 * the sumsl default is 100*machep
62 * absolute convergence
63 if (tolf.eq.0.0D0) tolf=1.0D-4
65 * relative convergence
66 if (rtolf.eq.0.0D0) rtolf=1.0D-4
68 * controls initial step size
70 * large vals of d correspond to small components of step
77 cd print *,'Calling SUMSL'
79 call x2xx(x,xx,nvar_restr)
80 call sumsl(nvar_restr,d,xx,func_restr,grad_restr,
81 & iv,liv,lv,v,idum,rdum,fdum)
84 call sumsl(nvar,d,x,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
88 cd print *,'Exit SUMSL; return code:',iretcode,' energy:',etot
89 write (iout,'(/a,i4,a,f15.5/)') 'SUMSL return code:',iv(1),
91 call var_to_geom(nvar,x)
93 c call etotal(energia(0))
95 c call enerprint(energia(0))
101 c----------------------------------------------------------------------------
102 subroutine ergastulum
103 implicit real*8 (a-h,o-z)
108 include 'COMMON.SETUP'
109 include 'COMMON.DERIV'
111 include 'COMMON.IOUNITS'
112 include 'COMMON.FFIELD'
113 include 'COMMON.INTERACT'
115 include 'COMMON.TIME1'
116 double precision z(maxres6),d_a_tmp(maxres6)
117 double precision edum(0:n_ene),time_order(0:10)
118 double precision Gcopy(maxres2,maxres2)
119 common /przechowalnia/ Gcopy
121 C Workers wait for variables and NF, and NFL from the boss
123 do while (iorder.ge.0)
124 c write (*,*) 'Processor',fg_rank,' CG group',kolor,
125 c & ' receives order from Master'
127 call MPI_Bcast(iorder,1,MPI_INTEGER,king,FG_COMM,IERR)
128 time_Bcast=time_Bcast+MPI_Wtime()-time00
129 if (icall.gt.4 .and. iorder.ge.0)
130 & time_order(iorder)=time_order(iorder)+MPI_Wtime()-time00
133 c & 'Processor',fg_rank,' completed receive MPI_BCAST order',iorder
134 if (iorder.eq.0) then
137 c write (2,*) "After etotal"
138 c write (2,*) "dimen",dimen," dimen3",dimen3
140 else if (iorder.eq.2) then
142 cmd call etotal_short(edum)
143 c write (2,*) "After etotal_short"
144 c write (2,*) "dimen",dimen," dimen3",dimen3
146 else if (iorder.eq.3) then
148 cmd call etotal_long(edum)
149 c write (2,*) "After etotal_long"
150 c write (2,*) "dimen",dimen," dimen3",dimen3
152 else if (iorder.eq.1) then
154 c write (2,*) "After sum_gradient"
155 c write (2,*) "dimen",dimen," dimen3",dimen3
157 else if (iorder.eq.4) then
158 cmd call ginv_mult(z,d_a_tmp)
159 else if (iorder.eq.5) then
160 c Setup MD things for a slave
161 dimen=(nct-nnt+1)+nside
162 dimen1=(nct-nnt)+(nct-nnt+1)
164 c write (2,*) "dimen",dimen," dimen3",dimen3
166 call int_bounds(dimen,igmult_start,igmult_end)
167 igmult_start=igmult_start-1
168 call MPI_Allgather(3*igmult_start,1,MPI_INTEGER,
169 & ng_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
170 my_ng_count=igmult_end-igmult_start
171 call MPI_Allgather(3*my_ng_count,1,MPI_INTEGER,ng_counts(0),1,
172 & MPI_INTEGER,FG_COMM,IERROR)
173 c write (2,*) "ng_start",(ng_start(i),i=0,nfgtasks-1)
174 c write (2,*) "ng_counts",(ng_counts(i),i=0,nfgtasks-1)
175 myginv_ng_count=maxres2*my_ng_count
176 c write (2,*) "igmult_start",igmult_start," igmult_end",
177 c & igmult_end," my_ng_count",my_ng_count
179 call MPI_Allgather(maxres2*igmult_start,1,MPI_INTEGER,
180 & nginv_start(0),1,MPI_INTEGER,FG_COMM,IERROR)
181 call MPI_Allgather(myginv_ng_count,1,MPI_INTEGER,
182 & nginv_counts(0),1,MPI_INTEGER,FG_COMM,IERROR)
183 c write (2,*) "nginv_start",(nginv_start(i),i=0,nfgtasks-1)
184 c write (2,*) "nginv_counts",(nginv_counts(i),i=0,nfgtasks-1)
186 c call MPI_Barrier(FG_COMM,IERROR)
188 call MPI_Scatterv(ginv(1,1),nginv_counts(0),
189 & nginv_start(0),MPI_DOUBLE_PRECISION,gcopy(1,1),
190 & myginv_ng_count,MPI_DOUBLE_PRECISION,king,FG_COMM,IERR)
192 time_scatter_ginv=time_scatter_ginv+MPI_Wtime()-time00
199 c write (2,*) "dimen",dimen," dimen3",dimen3
200 c write (2,*) "End MD setup"
202 c write (iout,*) "My chunk of ginv_block"
203 c call MATOUT2(my_ng_count,dimen3,maxres2,maxers2,ginv_block)
204 else if (iorder.eq.6) then
205 call int_from_cart1(.false.)
206 else if (iorder.eq.7) then
208 else if (iorder.eq.8) then
210 else if (iorder.eq.9) then
211 cmd call fricmat_mult(z,d_a_tmp)
212 else if (iorder.eq.10) then
213 cmd call setup_fricmat
216 write (*,*) 'Processor',fg_rank,' CG group',kolor,
217 & ' absolute rank',myrank,' leves ERGASTULUM.'
218 write(*,*)'Processor',fg_rank,' wait times for respective orders',
219 & (' order[',i,']',time_order(i),i=0,10)
223 ************************************************************************
224 subroutine func(n,x,nf,f,uiparm,urparm,ufparm)
225 implicit real*8 (a-h,o-z)
227 include 'COMMON.DERIV'
228 include 'COMMON.IOUNITS'
231 double precision energia(0:n_ene)
233 double precision ufparm
239 c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
243 cd print *,'func',nf,nfl,icg
244 call var_to_geom(n,x)
247 cd write (iout,*) 'ETOTAL called from FUNC'
248 call etotal(energia(0))
252 c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
253 c write (iout,*) 'f=',etot
258 ************************************************************************
259 subroutine func_restr(n,x,nf,f,uiparm,urparm,ufparm)
260 implicit real*8 (a-h,o-z)
262 include 'COMMON.DERIV'
263 include 'COMMON.IOUNITS'
266 double precision energia(0:n_ene)
268 double precision ufparm
274 c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
278 call var_to_geom_restr(n,x)
281 cd write (iout,*) 'ETOTAL called from FUNC'
282 call etotal(energia(0))
286 c write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
287 c write (iout,*) 'f=',etot
292 c-------------------------------------------------------
293 subroutine x2xx(x,xx,n)
294 implicit real*8 (a-h,o-z)
297 include 'COMMON.CHAIN'
298 include 'COMMON.INTERACT'
299 double precision xx(maxvar),x(maxvar)
309 if (mask_phi(i).eq.1) then
317 if (mask_theta(i).eq.1) then
325 if (itype(i).ne.10) then
327 if (mask_side(i).eq.1) then
340 subroutine xx2x(x,xx)
341 implicit real*8 (a-h,o-z)
344 include 'COMMON.CHAIN'
345 include 'COMMON.INTERACT'
346 double precision xx(maxvar),x(maxvar)
356 if (mask_phi(i).eq.1) then
364 if (mask_theta(i).eq.1) then
372 if (itype(i).ne.10) then
374 if (mask_side(i).eq.1) then
385 c----------------------------------------------------------
386 subroutine minim_dc(etot,iretcode,nfun)
387 implicit real*8 (a-h,o-z)
389 parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
393 include 'COMMON.SETUP'
394 include 'COMMON.IOUNITS'
397 include 'COMMON.MINIM'
398 include 'COMMON.CHAIN'
400 double precision minval,x(maxvar),d(maxvar),v(1:lv),xx(maxvar)
401 c common /przechowalnia/ v
403 double precision energia(0:n_ene)
404 external func_dc,grad_dc,fdum
405 logical not_done,change,reduce
406 double precision g(maxvar),f1
408 call deflt(2,iv,liv,lv,v)
409 * 12 means fresh start, dont call deflt
411 * max num of fun calls
412 if (maxfun.eq.0) maxfun=500
414 * max num of iterations
415 if (maxmin.eq.0) maxmin=1000
419 * selects output unit
420 if (print_min_ini+print_min_stat+print_min_res.gt.0) iv(21)=iout
421 * 1 means to print out result
423 * 1 means to print out summary stats
424 iv(23)=print_min_stat
425 * 1 means to print initial x and d
427 * min val for v(radfac) default is 0.1
429 * max val for v(radfac) default is 4.0
432 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
433 * the sumsl default is 0.1
435 * false conv if (act fnctn decrease) .lt. v(34)
436 * the sumsl default is 100*machep
438 * absolute convergence
439 if (tolf.eq.0.0D0) tolf=1.0D-4
441 * relative convergence
442 if (rtolf.eq.0.0D0) rtolf=1.0D-4
444 * controls initial step size
446 * large vals of d correspond to small components of step
459 if (ialph(i,1).gt.0) then
467 call sumsl(k,d,x,func_dc,grad_dc,iv,liv,lv,v,idum,rdum,fdum)
477 if (ialph(i,1).gt.0) then
488 cd call func_dc(k,x,nf,f,idum,rdum,fdum)
489 cd call grad_dc(k,x,nf,g,idum,rdum,fdum)
493 cd call func_dc(k,x,nf,f1,idum,rdum,fdum)
495 cd print '(i5,2f15.5)',i,g(i),(f1-f)/1.0D-5
504 subroutine func_dc(n,x,nf,f,uiparm,urparm,ufparm)
505 implicit real*8 (a-h,o-z)
510 include 'COMMON.SETUP'
511 include 'COMMON.DERIV'
512 include 'COMMON.IOUNITS'
514 include 'COMMON.CHAIN'
516 double precision energia(0:n_ene)
517 double precision ufparm
534 if (ialph(i,1).gt.0) then
544 call etotal(energia(0))
547 cd print *,'func_dc ',nf,nfl,f
552 subroutine grad_dc(n,x,nf,g,uiparm,urparm,ufparm)
553 implicit real*8 (a-h,o-z)
558 include 'COMMON.SETUP'
559 include 'COMMON.CHAIN'
560 include 'COMMON.DERIV'
562 include 'COMMON.INTERACT'
563 include 'COMMON.FFIELD'
565 include 'COMMON.IOUNITS'
568 double precision urparm(1)
569 dimension x(maxvar),g(maxvar)
575 cd print *,'grad_dc ',nf,nfl,nf-nfl+1,icg
576 if (nf-nfl+1) 20,30,40
577 20 call func_dc(n,x,nf,f,uiparm,urparm,ufparm)
591 if (ialph(i,1).gt.0) then
601 C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
613 if (ialph(i,1).gt.0) then