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