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