wham corrections
[unres.git] / source / wham / src-M / initialize_p.F.org
1       subroutine initialize
2
3 C Define constants and zero out tables.
4 C
5       implicit real*8 (a-h,o-z)
6       include 'DIMENSIONS'
7       include 'DIMENSIONS.ZSCOPT'
8 #ifdef MPI
9       include 'mpif.h'
10 #endif
11       include 'COMMON.IOUNITS'
12       include 'COMMON.CHAIN'
13       include 'COMMON.INTERACT'
14       include 'COMMON.GEO'
15       include 'COMMON.LOCAL'
16       include 'COMMON.TORSION'
17       include 'COMMON.FFIELD'
18       include 'COMMON.SBRIDGE'
19       include 'COMMON.MINIM' 
20       include 'COMMON.DERIV'
21       include "COMMON.WEIGHTS"
22       include "COMMON.NAMES"
23       include "COMMON.TIME1"
24 C
25 C The following is just to define auxiliary variables used in angle conversion
26 C
27       pi=4.0D0*datan(1.0D0)
28       dwapi=2.0D0*pi
29       dwapi3=dwapi/3.0D0
30       pipol=0.5D0*pi
31       deg2rad=pi/180.0D0
32       rad2deg=1.0D0/deg2rad
33       angmin=10.0D0*deg2rad
34 C
35 C Define I/O units.
36 C
37       inp=    1
38       iout=   2
39       ipdbin= 3
40       ipdb=   7
41       imol2=  4
42       igeom=  8
43       intin=  9
44       ithep= 11
45       irotam=12
46       itorp= 13
47       itordp= 23
48       ielep= 14
49       isidep=15 
50       iscpp=25
51       icbase=16
52       ifourier=20
53       istat= 17
54       ientin=18
55       ientout=19
56 C
57 C CSA I/O units (separated from others especially for Jooyoung)
58 C
59       icsa_rbank=30
60       icsa_seed=31
61       icsa_history=32
62       icsa_bank=33
63       icsa_bank1=34
64       icsa_alpha=35
65       icsa_alpha1=36
66       icsa_bankt=37
67       icsa_int=39
68       icsa_bank_reminimized=38
69       icsa_native_int=41
70       icsa_in=40
71 C
72 C Set default weights of the energy terms.
73 C
74       wlong=1.0D0
75       welec=1.0D0
76       wtor =1.0D0
77       wang =1.0D0
78       wscloc=1.0D0
79       wstrain=1.0D0
80 C
81 C Zero out tables.
82 C
83       ndih_constr=0
84       do i=1,maxres2
85         do j=1,3
86           c(j,i)=0.0D0
87           dc(j,i)=0.0D0
88         enddo
89       enddo
90       do i=1,maxres
91         do j=1,3
92           xloc(j,i)=0.0D0
93         enddo
94       enddo
95       do i=1,ntyp
96         do j=1,ntyp
97           aa(i,j)=0.0D0
98           bb(i,j)=0.0D0
99           augm(i,j)=0.0D0
100           sigma(i,j)=0.0D0
101           r0(i,j)=0.0D0
102           chi(i,j)=0.0D0
103         enddo
104         do j=1,2
105           bad(i,j)=0.0D0
106         enddo
107         chip(i)=0.0D0
108         alp(i)=0.0D0
109         sigma0(i)=0.0D0
110         sigii(i)=0.0D0
111         rr0(i)=0.0D0
112         a0thet(i)=0.0D0
113         do j=1,2
114           athet(j,i)=0.0D0
115           bthet(j,i)=0.0D0
116         enddo
117         do j=0,3
118           polthet(j,i)=0.0D0
119         enddo
120         do j=1,3
121           gthet(j,i)=0.0D0
122         enddo
123         theta0(i)=0.0D0
124         sig0(i)=0.0D0
125         sigc0(i)=0.0D0
126         do j=1,maxlob
127           bsc(j,i)=0.0D0
128           do k=1,3
129             censc(k,j,i)=0.0D0
130           enddo
131           do k=1,3
132             do l=1,3
133               gaussc(l,k,j,i)=0.0D0
134             enddo
135           enddo
136           nlob(i)=0
137         enddo
138       enddo
139       nlob(ntyp1)=0
140       dsc(ntyp1)=0.0D0
141       do i=1,maxtor
142         itortyp(i)=0
143         do j=1,maxtor
144           do k=1,maxterm
145             v1(k,j,i)=0.0D0
146             v2(k,j,i)=0.0D0
147           enddo
148         enddo
149       enddo
150       do i=1,maxres
151         itype(i)=0
152         itel(i)=0
153       enddo
154 C Initialize the bridge arrays
155       ns=0
156       nss=0 
157       nhpb=0
158       do i=1,maxss
159         iss(i)=0
160       enddo
161       do i=1,maxdim
162         dhpb(i)=0.0D0
163       enddo
164       do i=1,maxres
165         ihpb(i)=0
166         jhpb(i)=0
167       enddo
168 C
169 C Initialize timing.
170 C
171       call set_timers
172 C
173 C Initialize variables used in minimization.
174 C   
175 c     maxfun=5000
176 c     maxit=2000
177       maxfun=500
178       maxit=200
179       tolf=1.0D-2
180       rtolf=5.0D-4
181
182 C Initialize the variables responsible for the mode of gradient storage.
183 C
184       nfl=0
185       icg=1
186       do i=1,14
187         do j=1,14
188           if (print_order(i).eq.j) then
189             iw(print_order(i))=j
190             goto 1121
191           endif
192         enddo
193 1121    continue
194       enddo
195       calc_grad=.false.
196 C Set timers and counters for the respective routines
197       t_func = 0.0d0
198       t_grad = 0.0d0
199       t_fhel = 0.0d0
200       t_fbet = 0.0d0
201       t_ghel = 0.0d0
202       t_gbet = 0.0d0
203       t_viol = 0.0d0
204       t_gviol = 0.0d0
205       n_func = 0
206       n_grad = 0
207       n_fhel = 0
208       n_fbet = 0
209       n_ghel = 0
210       n_gbet = 0
211       n_viol = 0
212       n_gviol = 0
213       n_map = 0
214       return
215       end
216 c-------------------------------------------------------------------------
217       block data nazwy
218       implicit real*8 (a-h,o-z)
219       include 'DIMENSIONS'
220       include 'DIMENSIONS.ZSCOPT'
221       include 'COMMON.NAMES'
222       include 'COMMON.FFIELD'
223       data restyp /
224      &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
225      &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','D'/
226       data onelet /
227      &'C','M','F','I','L','V','W','Y','A','G','T',
228      &'S','Q','N','E','D','H','R','K','P','X'/
229       data potname /'LJ','LJK','BP','GB','GBV'/
230       data ename / 
231      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
232      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
233      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EVDW2_14",2*" "/
234       data wname /
235      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
236      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
237      &   "SCAL14",2*" "/
238 #ifdef SCP14
239       data nprint_ene /15/
240       data print_order /1,2,3,11,12,13,14,4,5,6,7,8,9,10,16,0/
241 #else
242       data nprint_ene /14/
243       data print_order /1,2,3,11,12,13,14,4,5,6,7,8,9,10,3*0/
244 #endif
245       end 
246 c---------------------------------------------------------------------------
247       subroutine init_int_table
248       implicit real*8 (a-h,o-z)
249       include 'DIMENSIONS'
250       include 'DIMENSIONS.ZSCOPT'
251 #ifdef MPI
252       include 'mpif.h'
253 #endif
254 #ifdef MP
255       include 'COMMON.INFO'
256 #endif
257       include 'COMMON.CHAIN'
258       include 'COMMON.INTERACT'
259       include 'COMMON.LOCAL'
260       include 'COMMON.SBRIDGE'
261       include 'COMMON.IOUNITS'
262       logical scheck,lprint
263 #ifdef MPL
264       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
265      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
266 C... Determine the numbers of start and end SC-SC interaction 
267 C... to deal with by current processor.
268       lprint=.false.
269       if (lprint)
270      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
271       n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
272       MyRank=MyID-(MyGroup-1)*fgProcs
273       call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
274       if (lprint)
275      &  write (iout,*) 'Processor',MyID,' MyRank',MyRank,
276      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
277      &  ' my_sc_inde',my_sc_inde
278       ind_sctint=0
279       iatsc_s=0
280       iatsc_e=0
281 #endif
282       lprint=.false.
283       do i=1,maxres
284         nint_gr(i)=0
285         nscp_gr(i)=0
286         do j=1,maxint_gr
287           istart(i,1)=0
288           iend(i,1)=0
289           ielstart(i)=0
290           ielend(i)=0
291           iscpstart(i,1)=0
292           iscpend(i,1)=0    
293         enddo
294       enddo
295       ind_scint=0
296       ind_scint_old=0
297 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
298 cd   &   (ihpb(i),jhpb(i),i=1,nss)
299       do i=nnt,nct-1
300         scheck=.false.
301         do ii=1,nss
302           if (ihpb(ii).eq.i+nres) then
303             scheck=.true.
304             jj=jhpb(ii)-nres
305             goto 10
306           endif
307         enddo
308    10   continue
309 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
310         if (scheck) then
311           if (jj.eq.i+1) then
312 #ifdef MPL
313             write (iout,*) 'jj=i+1'
314             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
315      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
316 #else
317             nint_gr(i)=1
318             istart(i,1)=i+2
319             iend(i,1)=nct
320 #endif
321           else if (jj.eq.nct) then
322 #ifdef MPL
323             write (iout,*) 'jj=nct'
324             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
325      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
326 #else
327             nint_gr(i)=1
328             istart(i,1)=i+1
329             iend(i,1)=nct-1
330 #endif
331           else
332 #ifdef MPL
333             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
334      & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
335             ii=nint_gr(i)+1
336             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
337      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
338 #else
339             nint_gr(i)=2
340             istart(i,1)=i+1
341             iend(i,1)=jj-1
342             istart(i,2)=jj+1
343             iend(i,2)=nct
344 #endif
345           endif
346         else
347 #ifdef MPL
348           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
349      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
350 #else
351           nint_gr(i)=1
352           istart(i,1)=i+1
353           iend(i,1)=nct
354           ind_scint=int_scint+nct-i
355 #endif
356         endif
357 #ifdef MPL
358         ind_scint_old=ind_scint
359 #endif
360       enddo
361    12 continue
362 #ifndef MPL
363       iatsc_s=nnt
364       iatsc_e=nct-1
365 #endif
366 #ifdef MPL
367       if (lprint) then
368         write (iout,*) 'Processor',MyID,' Group',MyGroup
369         write (iout,*) 'iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
370       endif
371 #endif
372       if (lprint) then
373       write (iout,'(a)') 'Interaction array:'
374       do i=iatsc_s,iatsc_e
375         write (iout,'(i3,2(2x,2i3))') 
376      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
377       enddo
378       endif
379       ispp=2
380 #ifdef MPL
381 C Now partition the electrostatic-interaction array
382       npept=nct-nnt
383       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
384       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
385       if (lprint)
386      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
387      &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
388      &               ' my_ele_inde',my_ele_inde
389       iatel_s=0
390       iatel_e=0
391       ind_eleint=0
392       ind_eleint_old=0
393       do i=nnt,nct-3
394         ijunk=0
395         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
396      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
397       enddo ! i 
398    13 continue
399 #else
400       iatel_s=nnt
401       iatel_e=nct-3
402       do i=iatel_s,iatel_e
403         ielstart(i)=i+2
404         ielend(i)=nct-1
405       enddo
406 #endif
407       if (lprint) then
408         write (iout,'(a)') 'Electrostatic interaction array:'
409         do i=iatel_s,iatel_e
410           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
411         enddo
412       endif ! lprint
413 c     iscp=3
414       iscp=2
415 C Partition the SC-p interaction array
416 #ifdef MPL
417       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
418       call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
419       if (lprint)
420      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
421      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
422      &               ' my_scp_inde',my_scp_inde
423       iatscp_s=0
424       iatscp_e=0
425       ind_scpint=0
426       ind_scpint_old=0
427       do i=nnt,nct-1
428         if (i.lt.nnt+iscp) then
429 cd        write (iout,*) 'i.le.nnt+iscp'
430           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
431      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
432      &      iscpend(i,1),*14)
433         else if (i.gt.nct-iscp) then
434 cd        write (iout,*) 'i.gt.nct-iscp'
435           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
436      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
437      &      iscpend(i,1),*14)
438         else
439           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
440      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
441      &      iscpend(i,1),*14)
442           ii=nscp_gr(i)+1
443           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
444      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
445      &      iscpend(i,ii),*14)
446         endif
447       enddo ! i
448    14 continue
449 #else
450       iatscp_s=nnt
451       iatscp_e=nct-1
452       do i=nnt,nct-1
453         if (i.lt.nnt+iscp) then
454           nscp_gr(i)=1
455           iscpstart(i,1)=i+iscp
456           iscpend(i,1)=nct
457         elseif (i.gt.nct-iscp) then
458           nscp_gr(i)=1
459           iscpstart(i,1)=nnt
460           iscpend(i,1)=i-iscp
461         else
462           nscp_gr(i)=2
463           iscpstart(i,1)=nnt
464           iscpend(i,1)=i-iscp
465           iscpstart(i,2)=i+iscp
466           iscpend(i,2)=nct
467         endif 
468       enddo ! i
469 #endif
470       if (lprint) then
471         write (iout,'(a)') 'SC-p interaction array:'
472         do i=iatscp_s,iatscp_e
473           write (iout,'(i3,2(2x,2i3))') 
474      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
475         enddo
476       endif ! lprint
477 C Partition local interactions
478 #ifdef MPL
479       call int_bounds(nres-2,loc_start,loc_end)
480       loc_start=loc_start+1
481       loc_end=loc_end+1
482       call int_bounds(nres-2,ithet_start,ithet_end)
483       ithet_start=ithet_start+2
484       ithet_end=ithet_end+2
485       call int_bounds(nct-nnt-2,iphi_start,iphi_end) 
486       iphi_start=iphi_start+nnt+2
487       iphi_end=iphi_end+nnt+2
488       if (lprint) then 
489         write (iout,*) 'Processor:',MyID,
490      & ' loc_start',loc_start,' loc_end',loc_end,
491      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
492      & ' iphi_start',iphi_start,' iphi_end',iphi_end
493         write (*,*) 'Processor:',MyID,
494      & ' loc_start',loc_start,' loc_end',loc_end,
495      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
496      & ' iphi_start',iphi_start,' iphi_end',iphi_end
497       endif
498       if (fgprocs.gt.1 .and. MyID.eq.BossID) then
499         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
500      & nele_int_tot,' electrostatic and ',nscp_int_tot,
501      & ' SC-p interactions','were distributed among',fgprocs,
502      & ' fine-grain processors.'
503       endif
504 #else
505       loc_start=2
506       loc_end=nres-1
507       ithet_start=3 
508       ithet_end=nres
509       iphi_start=nnt+3
510       iphi_end=nct
511 #endif
512       return
513       end 
514 c---------------------------------------------------------------------------
515       subroutine int_partition(int_index,lower_index,upper_index,atom,
516      & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
517       implicit real*8 (a-h,o-z)
518       include 'DIMENSIONS'
519       include 'COMMON.IOUNITS'
520       integer int_index,lower_index,upper_index,atom,at_start,at_end,
521      & first_atom,last_atom,int_gr,jat_start,jat_end
522       logical lprn
523       lprn=.false.
524       if (lprn) write (iout,*) 'int_index=',int_index
525       int_index_old=int_index
526       int_index=int_index+last_atom-first_atom+1
527       if (lprn) 
528      &   write (iout,*) 'int_index=',int_index,
529      &               ' int_index_old',int_index_old,
530      &               ' lower_index=',lower_index,
531      &               ' upper_index=',upper_index,
532      &               ' atom=',atom,' first_atom=',first_atom,
533      &               ' last_atom=',last_atom
534       if (int_index.ge.lower_index) then
535         int_gr=int_gr+1
536         if (at_start.eq.0) then
537           at_start=atom
538           jat_start=first_atom-1+lower_index-int_index_old
539         else
540           jat_start=first_atom
541         endif
542         if (lprn) write (iout,*) 'jat_start',jat_start
543         if (int_index.ge.upper_index) then
544           at_end=atom
545           jat_end=first_atom-1+upper_index-int_index_old
546           return1
547         else
548           jat_end=last_atom
549         endif
550         if (lprn) write (iout,*) 'jat_end',jat_end
551       endif
552       return
553       end
554 c------------------------------------------------------------------------------
555       subroutine hpb_partition
556       implicit real*8 (a-h,o-z)
557       include 'DIMENSIONS'
558       include 'COMMON.SBRIDGE'
559       include 'COMMON.IOUNITS'
560 #ifdef MPL
561       include 'COMMON.INFO'
562       call int_bounds(nhpb,link_start,link_end)
563 #else
564       link_start=1
565       link_end=nhpb
566 #endif
567 cd    write (iout,*) 'Processor',MyID,' MyRank',MyRank,
568 cd   &  ' nhpb',nhpb,' link_start=',link_start,
569 cd   &  ' link_end',link_end
570       return
571       end