6/11/2012 by Adam
[unres.git] / source / unres / src_CSA / minim_jlee.F
1       subroutine minim_jlee
2 c  controls minimization and sorting routines
3       implicit real*8 (a-h,o-z)
4       include 'DIMENSIONS'
5       parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2)) 
6       include 'COMMON.VAR'
7       include 'COMMON.IOUNITS'
8       include 'COMMON.MINIM'
9       include 'COMMON.CONTROL'
10       include 'mpif.h'
11       external func,gradient,fdum
12       real ran1,ran2,ran3
13       include 'COMMON.SETUP'
14       include 'COMMON.GEO'
15       include 'COMMON.FFIELD'
16       include 'COMMON.SBRIDGE'
17       include 'COMMON.DISTFIT'
18       include 'COMMON.CHAIN'
19       dimension muster(mpi_status_size)
20       dimension var(maxvar),erg(mxch*(mxch+1)/2+1)
21       dimension var2(maxvar)
22       integer iffr(maxres),ihpbt(maxdim),jhpbt(maxdim)
23       double precision d(maxvar),v(1:lv+1),garbage(maxvar)                     
24       double precision energia(0:n_ene),time0s,time1s
25       dimension indx(9),info(12)
26       dimension iv(liv)                                               
27       dimension idum(1),rdum(1)
28       dimension icont(2,maxcont)
29       logical check_var,fail
30       integer iloop(2)
31       common /przechowalnia/ v
32       data rad /1.745329252d-2/
33 c  receive # of start
34 !      print *,'Processor',me,' calling MINIM_JLEE maxfun',maxfun,
35 !     &   ' maxmin',maxmin,' tolf',tolf,' rtolf',rtolf
36       nhpb0=nhpb
37    10 continue
38       time0s=MPI_WTIME()
39 c      print *, 'MINIM_JLEE: ',me,' is waiting'
40       call mpi_recv(info,12,mpi_integer,king,idint,CG_COMM,
41      *              muster,ierr)
42       time1s=MPI_WTIME()
43       write (iout,'(a12,f10.4,a4)')'Waiting for ',time1s-time0s,' sec'
44       call flush(iout)
45       if (info(1).eq.0.and.info(2).eq.-2) then
46 cd       write (iout,*) 'Parallel tmscore for refresh bank'
47 cd       call flush(iout)
48        call refresh_bank_worker_tmscore(var)
49        goto 10
50       endif
51        n=info(1)
52 c      print *, 'MINIM_JLEE: ',me,' received: ',n
53       
54 crc      if (ierr.ne.0) go to 100
55 c  if # = 0, return
56       if (n.eq.0) then 
57         write (iout,*) 'Finishing minim_jlee - signal',n,' from master'
58         call flush(iout)
59         return
60       endif
61
62       nfun=0
63       IF (n.lt.0) THEN
64        call mpi_recv(var,nvar,mpi_double_precision,
65      *              king,idreal,CG_COMM,muster,ierr)
66        call mpi_recv(iffr,nres,mpi_integer,
67      *              king,idint,CG_COMM,muster,ierr)
68        call mpi_recv(var2,nvar,mpi_double_precision,
69      *              king,idreal,CG_COMM,muster,ierr)
70       ELSE
71 c  receive initial values of variables
72        call mpi_recv(var,nvar,mpi_double_precision,
73      *              king,idreal,CG_COMM,muster,ierr)
74 crc       if (ierr.ne.0) go to 100
75       ENDIF
76
77       if(vdisulf.and.info(2).ne.-1) then
78         if(info(4).ne.0)then
79            call mpi_recv(ihpbt,info(4),mpi_integer,
80      *              king,idint,CG_COMM,muster,ierr)
81            call mpi_recv(jhpbt,info(4),mpi_integer,
82      *              king,idint,CG_COMM,muster,ierr)
83         endif
84       endif  
85
86       IF (n.lt.0) THEN
87         n=-n     
88         nhpb=nhpb0
89         link_start=1
90         link_end=nhpb
91         call init_int_table
92         call contact_cp(var,var2,iffr,nfun,n)
93       ENDIF
94
95       if(vdisulf.and.info(2).ne.-1) then
96         nss=0
97         if(info(4).ne.0)then
98 cd          write(iout,*) 'SS=',info(4),'N=',info(1),'IT=',info(2)
99           call var_to_geom(nvar,var)
100           call chainbuild
101           do i=1,info(4)
102            if (dist(ihpbt(i),jhpbt(i)).lt.7.0) then
103             nss=nss+1
104             ihpb(nss)=ihpbt(i)
105             jhpb(nss)=jhpbt(i)
106 cd            write(iout,*) 'SS  mv=',info(3),
107 cd     &         ihpb(nss)-nres,jhpb(nss)-nres,
108 cd     &         dist(ihpb(nss),jhpb(nss))
109             dhpb(nss)=dbr
110             forcon(nss)=fbr
111            else
112 cd            write(iout,*) 'rm SS  mv=',info(3),
113 cd     &         ihpbt(i)-nres,jhpbt(i)-nres,dist(ihpbt(i),jhpbt(i))
114            endif
115           enddo
116         endif
117         nhpb=nss
118         link_start=1
119         link_end=nhpb
120         call init_int_table
121       endif
122
123       if (info(3).eq.14) then
124        write(iout,*) 'calling local_move',info(7),info(8)
125        call local_move_init(.false.)
126        call var_to_geom(nvar,var)
127        call local_move(info(7),info(8),20d0,50d0)
128        call geom_to_var(nvar,var)
129       endif
130
131
132       if (info(3).eq.16) then
133        write(iout,*) 'calling beta_slide',info(7),info(8),
134      &            info(10), info(11), info(12)
135        call var_to_geom(nvar,var)
136        call beta_slide(info(7),info(8),info(10),info(11),info(12)
137      &                                                     ,nfun,n)
138        call geom_to_var(nvar,var)
139       endif
140
141
142       if (info(3).eq.17) then
143        write(iout,*) 'calling beta_zip',info(7),info(8)
144        call var_to_geom(nvar,var)
145        call beta_zip(info(7),info(8),nfun,n)
146        call geom_to_var(nvar,var)
147       endif
148
149
150 crc overlap test
151
152       if (overlapsc) then 
153
154           call var_to_geom(nvar,var)
155           call chainbuild
156           call etotal(energia(0))
157           nfun=nfun+1
158           if (energia(1).eq.1.0d20) then  
159             info(3)=-info(3) 
160             write (iout,'(a,1pe14.5)')'#OVERLAP evdw=1d20',energia(1)
161             call overlap_sc(fail)
162             if(.not.fail) then
163              call geom_to_var(nvar,var)
164              call etotal(energia(0))
165              nfun=nfun+1
166              write (iout,'(a,1pe14.5)')'#OVERLAP evdw after',energia(1)
167             else
168              v(10)=1.0d20
169              iv(1)=-1
170              goto 201
171             endif
172           endif
173       endif 
174
175       if (searchsc) then 
176           call var_to_geom(nvar,var)
177           call sc_move(2,nres-1,1,10d0,nft_sc,etot)
178           call geom_to_var(nvar,var)
179 cd          write(iout,*) 'sc_move',nft_sc,etot
180       endif 
181
182       if (check_var(var,info)) then 
183            v(10)=1.0d21
184            iv(1)=6
185            goto 201
186       endif
187
188
189 crc        
190
191 !      write (iout,*) 'MINIM_JLEE: Processor',me,' nvar',nvar
192 !      write (iout,'(8f10.4)') (var(i),i=1,nvar)
193 !      write (*,*) 'MINIM_JLEE: Processor',me,' received nvar',nvar
194 !      write (*,'(8f10.4)') (var(i),i=1,nvar)
195
196        do i=1,nvar
197          garbage(i)=var(i)
198        enddo
199
200       call deflt(2,iv,liv,lv,v)                                         
201 * 12 means fresh start, dont call deflt                                 
202       iv(1)=12                                                          
203 * max num of fun calls                                                  
204       if (maxfun.eq.0) maxfun=500
205       iv(17)=maxfun
206 * max num of iterations                                                 
207       if (maxmin.eq.0) maxmin=1000
208       iv(18)=maxmin
209 * controls output                                                       
210       iv(19)=2                                                          
211 * selects output unit                                                   
212 cd      iv(21)=iout                                                       
213       iv(21)=0
214 * 1 means to print out result                                           
215       iv(22)=0                                                          
216 cd        iv(22)=1
217 * 1 means to print out summary stats                                    
218       iv(23)=0                                                          
219 * 1 means to print initial x and d                                      
220       iv(24)=0                                                          
221
222 c      if(me.eq.3.and.n.eq.255) then
223 c       print *,' CHUJ: stoi'
224 c       iv(21)=6
225 c       iv(22)=1
226 c       iv(23)=1
227 c       iv(24)=1                                                          
228 c      endif
229
230 * min val for v(radfac) default is 0.1                                  
231       v(24)=0.1D0                                                       
232 * max val for v(radfac) default is 4.0                                  
233       v(25)=2.0D0                                                       
234 c     v(25)=4.0D0                                                       
235 * check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)    
236 * the sumsl default is 0.1                                              
237       v(26)=0.1D0
238 * false conv if (act fnctn decrease) .lt. v(34)                         
239 * the sumsl default is 100*machep                                       
240       v(34)=v(34)/100.0D0                                               
241 * absolute convergence                                                  
242       if (tolf.eq.0.0D0) tolf=1.0D-4
243       v(31)=tolf
244 * relative convergence                                                  
245       if (rtolf.eq.0.0D0) rtolf=1.0D-4
246       v(32)=rtolf
247 * controls initial step size                                            
248        v(35)=1.0D-1                                                    
249 * large vals of d correspond to small components of step                
250       do i=1,nphi
251         d(i)=1.0D-1
252       enddo
253       do i=nphi+1,nvar
254         d(i)=1.0D-1
255       enddo
256 c  minimize energy
257 !      write (iout,*) 'Processor',me,' nvar',nvar
258 !      write (iout,*) 'Variables BEFORE minimization:'
259 !      write (iout,'(8f10.4)') (rad2deg*var(i),i=1,nvar)
260
261 c      print *, 'MINIM_JLEE: ',me,' before SUMSL '
262
263       call func(nvar,var,nf,eee,idum,rdum,fdum)
264       nfun=nfun+1
265       if(eee.ge.1.0d20) then
266 c       print *,'MINIM_JLEE: ',me,' CHUJ NASTAPIL'
267 c       print *,' energy before SUMSL =',eee
268 c       print *,' aborting local minimization'
269        iv(1)=-1
270        v(10)=eee
271        go to 201
272       endif
273
274 ct      time0s=MPI_WTIME()
275       call sumsl(nvar,d,var,func,gradient,iv,liv,lv,v,idum,rdum,fdum)
276 ct      write(iout,*) 'sumsl time=',MPI_WTIME()-time0s,iv(7),v(10)
277 c      print *, 'MINIM_JLEE: ',me,' after SUMSL '
278
279 c  find which conformation was returned from sumsl
280         nfun=nfun+iv(7)
281 !      print *,'Processor',me,' iv(17)',iv(17),' iv(18)',iv(18),' nf',nf,
282 !     & ' retcode',iv(1),' energy',v(10),' tolf',v(31),' rtolf',v(32)
283 c        if (iv(1).ne.4 .or. nf.le.1) then
284 c          write (*,*) 'Processor',me,' something bad in SUMSL',iv(1),nf
285 c         write (*,*) 'Initial Variables'
286 c         write (*,'(8f10.4)') (rad2deg*garbage(i),i=1,nvar)
287 c         write (*,*) 'Variables'
288 c         write (*,'(8f10.4)') (rad2deg*var(i),i=1,nvar)
289 c         write (*,*) 'Vector d'
290 c         write (*,'(8f10.4)') (d(i),i=1,nvar)
291 c         write (iout,*) 'Processor',me,' something bad in SUMSL',
292 c    &       iv(1),nf
293 c         write (iout,*) 'Initial Variables'
294 c         write (iout,'(8f10.4)') (rad2deg*garbage(i),i=1,nvar)
295 c         write (iout,*) 'Variables'
296 c         write (iout,'(8f10.4)') (rad2deg*var(i),i=1,nvar)
297 c         write (iout,*) 'Vector d'
298 c         write (iout,'(8f10.4)') (d(i),i=1,nvar)
299 c        endif
300 c        if (nf.lt.iv(6)-1) then
301 c  recalculate intra- and interchain energies
302 c         call func(nvar,var,nf,v(10),iv,v,fdum)
303 c        else if (nf.eq.iv(6)-1) then
304 c  regenerate conformation
305 c         call var_to_geom(nvar,var)
306 c         call chainbuild
307 c        endif
308 c  change origin and axes to standard ECEPP format
309 c        call var_to_geom(nvar,var)
310 !      write (iout,*) 'MINIM_JLEE after minim: Processor',me,' nvar',nvar
311 !      write (iout,'(8f10.4)') (var(i),i=1,nvar)
312 !      write (iout,*) 'Energy:',v(10)
313 c  send back output
314 c       print *, 'MINIM_JLEE: ',me,' minimized: ',n
315   201  continue
316         indx(1)=n
317 c return code: 6-gradient 9-number of ftn evaluation, etc
318         indx(2)=iv(1)
319 c total # of ftn evaluations (for iwf=0, it includes all minimizations).
320         indx(3)=nfun
321         indx(4)=info(2)
322         indx(5)=info(3)
323         indx(6)=nss
324         indx(7)=info(5)
325         indx(8)=info(6)
326         indx(9)=info(9)
327         call mpi_send(indx,9,mpi_integer,king,idint,CG_COMM,
328      *                 ierr)
329 c  send back energies
330 c al & cc
331 c calculate contact order
332 #ifdef CO_BIAS
333         call contact(.false.,ncont,icont,co)
334         erg(1)=v(10)-1.0d2*co
335 #else
336         erg(1)=v(10)
337 #endif
338         j=1
339         call mpi_send(erg,j,mpi_double_precision,king,idreal,
340      *                 CG_COMM,ierr)
341 #ifdef CO_BIAS
342         call mpi_send(co,j,mpi_double_precision,king,idreal,
343      *                 CG_COMM,ierr)
344 #endif
345 c  send back values of variables
346         call mpi_send(var,nvar,mpi_double_precision,
347      *               king,idreal,CG_COMM,ierr)
348 !        print * , 'MINIM_JLEE: Processor',me,' send erg and var '
349
350         if(vdisulf.and.info(2).ne.-1.and.nss.ne.0) then
351 cd          call intout
352 cd          call chainbuild
353 cd          call etotal(energia(0))
354 cd          etot=energia(0)
355 cd          call enerprint(energia(0))
356          call mpi_send(ihpb,nss,mpi_integer,
357      *               king,idint,CG_COMM,ierr)        
358          call mpi_send(jhpb,nss,mpi_integer,
359      *               king,idint,CG_COMM,ierr)        
360         endif
361
362         go to 10
363   100 print *, ' error in receiving message from emperor', me
364       call mpi_abort(mpi_comm_world,ierror,ierrcode)
365       return
366   200 print *, ' error in sending message to emperor'
367       call mpi_abort(mpi_comm_world,ierror,ierrcode)
368       return
369   300 print *, ' error in communicating with emperor'
370       call mpi_abort(mpi_comm_world,ierror,ierrcode)
371       return
372   956 format (' initial energy could not be calculated',41x)
373   957 format (80x)
374   965 format (' convergence code ',i2,' # of function calls ',
375      *  i4,' # of gradient calls ',i4,10x)
376   975 format (' energy ',1p,e12.4,' scaled gradient ',e11.3,32x)
377       end
378
379       logical function check_var(var,info)
380       implicit real*8 (a-h,o-z)
381       include 'DIMENSIONS'
382       include 'COMMON.VAR'
383       include 'COMMON.IOUNITS'
384       include 'COMMON.GEO'
385       include 'COMMON.SETUP'
386       dimension var(maxvar)
387       dimension info(3)
388 C AL -------
389        check_var=.false.
390        do i=nphi+ntheta+1,nphi+ntheta+nside
391 ! Check the side chain "valence" angles alpha
392          if (var(i).lt.1.0d-7) then 
393            write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
394            write (iout,*) 'Processor',me,'received bad variables!!!!'
395            write (iout,*) 'Variables'
396            write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
397            write (iout,*) 'Continuing calculations at this point',
398      & ' could destroy the results obtained so far... ABORTING!!!!!!'
399            write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)') 
400      &     'valence angle alpha',i-nphi-ntheta,var(i),
401      &     'n it',info(1),info(2),'mv ',info(3)
402            write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
403            write (*,*) 'Processor',me,'received bad variables!!!!'
404            write (*,*) 'Variables'
405            write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
406            write (*,*) 'Continuing calculations at this point',
407      & ' could destroy the results obtained so far... ABORTING!!!!!!'
408            write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)') 
409      &     'valence angle alpha',i-nphi-ntheta,var(i),
410      &     'n it',info(1),info(2),'mv ',info(3)
411            check_var=.true.
412            return
413          endif
414        enddo
415 ! Check the backbone "valence" angles theta
416        do i=nphi+1,nphi+ntheta
417          if (var(i).lt.1.0d-7) then
418            write (iout,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
419            write (iout,*) 'Processor',me,'received bad variables!!!!'
420            write (iout,*) 'Variables'
421            write (iout,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
422            write (iout,*) 'Continuing calculations at this point',
423      & ' could destroy the results obtained so far... ABORTING!!!!!!'
424            write (iout,'(a19,i5,f10.4,a4,2i4,a3,i3)') 
425      &     'valence angle theta',i-nphi,var(i),
426      &     'n it',info(1),info(2),'mv ',info(3)
427            write (*,*) 'CHUJ NASTAPIL ABSOLUTNY!!!!!!!!!!!!'
428            write (*,*) 'Processor',me,'received bad variables!!!!'
429            write (*,*) 'Variables'
430            write (*,'(8f10.4)') (rad2deg*var(j),j=1,nvar)
431            write (*,*) 'Continuing calculations at this point',
432      & ' could destroy the results obtained so far... ABORTING!!!!!!'
433            write (*,'(a19,i5,f10.4,a4,2i4,a3,i3)') 
434      &     'valence angle theta',i-nphi,var(i),
435      &     'n it',info(1),info(2),'mv ',info(3)
436            check_var=.true.
437            return
438          endif
439        enddo
440        return
441        end