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