fix 80 column in initialize_p.F
[unres.git] / source / unres / src_MD-M / initialize_p.F
1       block data
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4       include 'COMMON.MCM'
5       include 'COMMON.MD'
6       data MovTypID
7      &  /'pool','chain regrow','multi-bond','phi','theta','side chain',
8      &   'total'/
9 c Conversion from poises to molecular unit and the gas constant
10       data cPoise /2.9361d0/, Rb /0.001986d0/
11       end
12 c--------------------------------------------------------------------------
13       subroutine initialize
14
15 C Define constants and zero out tables.
16 C
17       implicit real*8 (a-h,o-z)
18       include 'DIMENSIONS'
19 #ifdef MPI
20       include 'mpif.h'
21 #endif
22 #ifndef ISNAN
23       external proc_proc
24 #ifdef WINPGI
25 cMS$ATTRIBUTES C ::  proc_proc
26 #endif
27 #endif
28       include 'COMMON.IOUNITS'
29       include 'COMMON.CHAIN'
30       include 'COMMON.INTERACT'
31       include 'COMMON.GEO'
32       include 'COMMON.LOCAL'
33       include 'COMMON.TORSION'
34       include 'COMMON.FFIELD'
35       include 'COMMON.SBRIDGE'
36       include 'COMMON.MCM'
37       include 'COMMON.MINIM' 
38       include 'COMMON.DERIV'
39       include 'COMMON.SPLITELE'
40 c Common blocks from the diagonalization routines
41       COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(400)
42       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
43       logical mask_r
44 c      real*8 text1 /'initial_i'/
45
46       mask_r=.false.
47 #ifndef ISNAN
48 c NaNQ initialization
49       i=-1
50       rr=dacos(100.0d0)
51 #ifdef WINPGI
52       idumm=proc_proc(rr,i)
53 #else
54       call proc_proc(rr,i)
55 #endif
56 #endif
57
58       kdiag=0
59       icorfl=0
60       iw=2
61 C
62 C The following is just to define auxiliary variables used in angle conversion
63 C
64       pi=4.0D0*datan(1.0D0)
65       dwapi=2.0D0*pi
66       dwapi3=dwapi/3.0D0
67       pipol=0.5D0*pi
68       deg2rad=pi/180.0D0
69       rad2deg=1.0D0/deg2rad
70       angmin=10.0D0*deg2rad
71 C
72 C Define I/O units.
73 C
74       inp=    1
75       iout=   2
76       ipdbin= 3
77       ipdb=   7
78       icart = 30
79       imol2=  4
80       igeom=  8
81       intin=  9
82       ithep= 11
83       irotam=12
84       itorp= 13
85       itordp= 23
86       ielep= 14
87       isidep=15 
88       iscpp=25
89       icbase=16
90       ifourier=20
91       istat= 17
92       irest1=55
93       irest2=56
94       iifrag=57
95       ientin=18
96       ientout=19
97       ibond = 28
98       isccor = 29
99 crc for write_rmsbank1  
100       izs1=21
101 cdr  include secondary structure prediction bias
102       isecpred=27
103 C
104 C CSA I/O units (separated from others especially for Jooyoung)
105 C
106       icsa_rbank=30
107       icsa_seed=31
108       icsa_history=32
109       icsa_bank=33
110       icsa_bank1=34
111       icsa_alpha=35
112       icsa_alpha1=36
113       icsa_bankt=37
114       icsa_int=39
115       icsa_bank_reminimized=38
116       icsa_native_int=41
117       icsa_in=40
118 crc for ifc error 118
119       icsa_pdb=42
120 C
121 C Set default weights of the energy terms.
122 C
123       wlong=1.0D0
124       welec=1.0D0
125       wtor =1.0D0
126       wang =1.0D0
127       wscloc=1.0D0
128       wstrain=1.0D0
129 C
130 C Zero out tables.
131 C
132 c      print '(a,$)','Inside initialize'
133 c      call memmon_print_usage()
134       do i=1,maxres2
135         do j=1,3
136           c(j,i)=0.0D0
137           dc(j,i)=0.0D0
138         enddo
139       enddo
140       do i=1,maxres
141         do j=1,3
142           xloc(j,i)=0.0D0
143         enddo
144       enddo
145       do i=1,ntyp
146         do j=1,ntyp
147           aa(i,j)=0.0D0
148           bb(i,j)=0.0D0
149           augm(i,j)=0.0D0
150           sigma(i,j)=0.0D0
151           r0(i,j)=0.0D0
152           chi(i,j)=0.0D0
153         enddo
154         do j=1,2
155           bad(i,j)=0.0D0
156         enddo
157         chip(i)=0.0D0
158         alp(i)=0.0D0
159         sigma0(i)=0.0D0
160         sigii(i)=0.0D0
161         rr0(i)=0.0D0
162         a0thet(i)=0.0D0
163         do j=1,2
164           athet(j,i)=0.0D0
165           bthet(j,i)=0.0D0
166         enddo
167         do j=0,3
168           polthet(j,i)=0.0D0
169         enddo
170         do j=1,3
171           gthet(j,i)=0.0D0
172         enddo
173         theta0(i)=0.0D0
174         sig0(i)=0.0D0
175         sigc0(i)=0.0D0
176         do j=1,maxlob
177           bsc(j,i)=0.0D0
178           do k=1,3
179             censc(k,j,i)=0.0D0
180           enddo
181           do k=1,3
182             do l=1,3
183               gaussc(l,k,j,i)=0.0D0
184             enddo
185           enddo
186           nlob(i)=0
187         enddo
188       enddo
189       nlob(ntyp1)=0
190       dsc(ntyp1)=0.0D0
191       do i=1,maxtor
192         itortyp(i)=0
193         do j=1,maxtor
194           do k=1,maxterm
195             v1(k,j,i)=0.0D0
196             v2(k,j,i)=0.0D0
197           enddo
198         enddo
199       enddo
200       do i=1,maxres
201         itype(i)=0
202         itel(i)=0
203       enddo
204 C Initialize the bridge arrays
205       ns=0
206       nss=0 
207       nhpb=0
208       do i=1,maxss
209         iss(i)=0
210       enddo
211       do i=1,maxdim
212         dhpb(i)=0.0D0
213       enddo
214       do i=1,maxres
215         ihpb(i)=0
216         jhpb(i)=0
217       enddo
218 C
219 C Initialize timing.
220 C
221       call set_timers
222 C
223 C Initialize variables used in minimization.
224 C   
225 c     maxfun=5000
226 c     maxit=2000
227       maxfun=500
228       maxit=200
229       tolf=1.0D-2
230       rtolf=5.0D-4
231
232 C Initialize the variables responsible for the mode of gradient storage.
233 C
234       nfl=0
235       icg=1
236 C
237 C Initialize constants used to split the energy into long- and short-range
238 C components
239 C
240       r_cut=2.0d0
241       rlamb=0.3d0
242 #ifndef SPLITELE
243       nprint_ene=nprint_ene-1
244 #endif
245       return
246       end
247 c-------------------------------------------------------------------------
248       block data nazwy
249       implicit real*8 (a-h,o-z)
250       include 'DIMENSIONS'
251       include 'COMMON.NAMES'
252       include 'COMMON.FFIELD'
253       data restyp /
254      &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
255      &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','D'/
256       data onelet /
257      &'C','M','F','I','L','V','W','Y','A','G','T',
258      &'S','Q','N','E','D','H','R','K','P','X'/
259       data potname /'LJ','LJK','BP','GB','GBV'/
260       data ename /
261      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
262      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
263      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB ","EVDWPP ",
264      &   "ESTR ","EVDW2_14 ","UCONST ", "      ","ESCCOR"/
265       data wname /
266      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
267      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
268      &   "WSTRAIN","WVDWPP","WBOND","SCAL14","     ","    ","WSCCOR"/
269       data nprint_ene /20/
270       data print_order/1,2,3,11,12,13,14,4,5,6,7,8,9,10,19,18,15,17,16,
271      & 21,0/
272       end 
273 c---------------------------------------------------------------------------
274       subroutine init_int_table
275       implicit real*8 (a-h,o-z)
276       include 'DIMENSIONS'
277 #ifdef MPI
278       include 'mpif.h'
279       integer blocklengths(15),displs(15)
280 #endif
281       include 'COMMON.CONTROL'
282       include 'COMMON.SETUP'
283       include 'COMMON.CHAIN'
284       include 'COMMON.INTERACT'
285       include 'COMMON.LOCAL'
286       include 'COMMON.SBRIDGE'
287       include 'COMMON.TORCNSTR'
288       include 'COMMON.IOUNITS'
289       include 'COMMON.DERIV'
290       include 'COMMON.CONTACTS'
291       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
292      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
293      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
294      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
295      & ielend_all(maxres,0:max_fg_procs-1),
296      & ntask_cont_from_all(0:max_fg_procs-1),
297      & itask_cont_from_all(0:max_fg_procs-1,0:max_fg_procs-1),
298      & ntask_cont_to_all(0:max_fg_procs-1),
299      & itask_cont_to_all(0:max_fg_procs-1,0:max_fg_procs-1)
300       integer FG_GROUP,CONT_FROM_GROUP,CONT_TO_GROUP
301       logical scheck,lprint,flag
302 #ifdef MPI
303       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
304      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
305 C... Determine the numbers of start and end SC-SC interaction 
306 C... to deal with by current processor.
307       do i=0,nfgtasks-1
308         itask_cont_from(i)=fg_rank
309         itask_cont_to(i)=fg_rank
310       enddo
311       lprint=.false.
312       if (lprint)
313      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
314       n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
315       call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
316       if (lprint)
317      &  write (iout,*) 'Processor',fg_rank,' CG group',kolor,
318      &  ' absolute rank',MyRank,
319      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
320      &  ' my_sc_inde',my_sc_inde
321       ind_sctint=0
322       iatsc_s=0
323       iatsc_e=0
324 #endif
325 c      lprint=.false.
326       do i=1,maxres
327         nint_gr(i)=0
328         nscp_gr(i)=0
329         do j=1,maxint_gr
330           istart(i,1)=0
331           iend(i,1)=0
332           ielstart(i)=0
333           ielend(i)=0
334           iscpstart(i,1)=0
335           iscpend(i,1)=0    
336         enddo
337       enddo
338       ind_scint=0
339       ind_scint_old=0
340 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
341 cd   &   (ihpb(i),jhpb(i),i=1,nss)
342       do i=nnt,nct-1
343         scheck=.false.
344         do ii=1,nss
345           if (ihpb(ii).eq.i+nres) then
346             scheck=.true.
347             jj=jhpb(ii)-nres
348             goto 10
349           endif
350         enddo
351    10   continue
352 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
353         if (scheck) then
354           if (jj.eq.i+1) then
355 #ifdef MPI
356 c            write (iout,*) 'jj=i+1'
357             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
358      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
359 #else
360             nint_gr(i)=1
361             istart(i,1)=i+2
362             iend(i,1)=nct
363 #endif
364           else if (jj.eq.nct) then
365 #ifdef MPI
366 c            write (iout,*) 'jj=nct'
367             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
368      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
369 #else
370             nint_gr(i)=1
371             istart(i,1)=i+1
372             iend(i,1)=nct-1
373 #endif
374           else
375 #ifdef MPI
376             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
377      & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
378             ii=nint_gr(i)+1
379             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
380      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
381 #else
382             nint_gr(i)=2
383             istart(i,1)=i+1
384             iend(i,1)=jj-1
385             istart(i,2)=jj+1
386             iend(i,2)=nct
387 #endif
388           endif
389         else
390 #ifdef MPI
391           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
392      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
393 #else
394           nint_gr(i)=1
395           istart(i,1)=i+1
396           iend(i,1)=nct
397           ind_scint=ind_scint+nct-i
398 #endif
399         endif
400 #ifdef MPI
401         ind_scint_old=ind_scint
402 #endif
403       enddo
404    12 continue
405 #ifndef MPI
406       iatsc_s=nnt
407       iatsc_e=nct-1
408 #endif
409 #ifdef MPI
410       if (lprint) write (*,*) 'Processor',fg_rank,' CG Group',kolor,
411      &   ' absolute rank',myrank,' iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
412 #endif
413       if (lprint) then
414       write (iout,'(a)') 'Interaction array:'
415       do i=iatsc_s,iatsc_e
416         write (iout,'(i3,2(2x,2i3))') 
417      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
418       enddo
419       endif
420       ispp=4
421 #ifdef MPI
422 C Now partition the electrostatic-interaction array
423       npept=nct-nnt
424       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
425       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
426       if (lprint)
427      & write (*,*) 'Processor',fg_rank,' CG group',kolor,
428      &  ' absolute rank',MyRank,
429      &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
430      &               ' my_ele_inde',my_ele_inde
431       iatel_s=0
432       iatel_e=0
433       ind_eleint=0
434       ind_eleint_old=0
435       do i=nnt,nct-3
436         ijunk=0
437         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
438      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
439       enddo ! i 
440    13 continue
441       if (iatel_s.eq.0) iatel_s=1
442       nele_int_tot_vdw=(npept-2)*(npept-2+1)/2
443 c      write (iout,*) "nele_int_tot_vdw",nele_int_tot_vdw
444       call int_bounds(nele_int_tot_vdw,my_ele_inds_vdw,my_ele_inde_vdw)
445 c      write (iout,*) "my_ele_inds_vdw",my_ele_inds_vdw,
446 c     & " my_ele_inde_vdw",my_ele_inde_vdw
447       ind_eleint_vdw=0
448       ind_eleint_vdw_old=0
449       iatel_s_vdw=0
450       iatel_e_vdw=0
451       do i=nnt,nct-3
452         ijunk=0
453         call int_partition(ind_eleint_vdw,my_ele_inds_vdw,
454      &    my_ele_inde_vdw,i,
455      &    iatel_s_vdw,iatel_e_vdw,i+2,nct-1,ijunk,ielstart_vdw(i),
456      &    ielend_vdw(i),*15)
457 c        write (iout,*) i," ielstart_vdw",ielstart_vdw(i),
458 c     &   " ielend_vdw",ielend_vdw(i)
459       enddo ! i 
460       if (iatel_s_vdw.eq.0) iatel_s_vdw=1
461    15 continue
462 #else
463       iatel_s=nnt
464       iatel_e=nct-5
465       do i=iatel_s,iatel_e
466         ielstart(i)=i+4
467         ielend(i)=nct-1
468       enddo
469       iatel_s_vdw=nnt
470       iatel_e_vdw=nct-3
471       do i=iatel_s_vdw,iatel_e_vdw
472         ielstart_vdw(i)=i+2
473         ielend_vdw(i)=nct-1
474       enddo
475 #endif
476       if (lprint) then
477         write (*,'(a)') 'Processor',fg_rank,' CG group',kolor,
478      &  ' absolute rank',MyRank
479         write (iout,*) 'Electrostatic interaction array:'
480         do i=iatel_s,iatel_e
481           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
482         enddo
483       endif ! lprint
484 c     iscp=3
485       iscp=2
486 C Partition the SC-p interaction array
487 #ifdef MPI
488       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
489       call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
490       if (lprint) write (iout,*) 'Processor',fg_rank,' CG group',kolor,
491      &  ' absolute rank',myrank,
492      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
493      &               ' my_scp_inde',my_scp_inde
494       iatscp_s=0
495       iatscp_e=0
496       ind_scpint=0
497       ind_scpint_old=0
498       do i=nnt,nct-1
499         if (i.lt.nnt+iscp) then
500 cd        write (iout,*) 'i.le.nnt+iscp'
501           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
502      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
503      &      iscpend(i,1),*14)
504         else if (i.gt.nct-iscp) then
505 cd        write (iout,*) 'i.gt.nct-iscp'
506           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
507      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
508      &      iscpend(i,1),*14)
509         else
510           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
511      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
512      &      iscpend(i,1),*14)
513           ii=nscp_gr(i)+1
514           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
515      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
516      &      iscpend(i,ii),*14)
517         endif
518       enddo ! i
519    14 continue
520 #else
521       iatscp_s=nnt
522       iatscp_e=nct-1
523       do i=nnt,nct-1
524         if (i.lt.nnt+iscp) then
525           nscp_gr(i)=1
526           iscpstart(i,1)=i+iscp
527           iscpend(i,1)=nct
528         elseif (i.gt.nct-iscp) then
529           nscp_gr(i)=1
530           iscpstart(i,1)=nnt
531           iscpend(i,1)=i-iscp
532         else
533           nscp_gr(i)=2
534           iscpstart(i,1)=nnt
535           iscpend(i,1)=i-iscp
536           iscpstart(i,2)=i+iscp
537           iscpend(i,2)=nct
538         endif 
539       enddo ! i
540 #endif
541       if (lprint) then
542         write (iout,'(a)') 'SC-p interaction array:'
543         do i=iatscp_s,iatscp_e
544           write (iout,'(i3,2(2x,2i3))') 
545      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
546         enddo
547       endif ! lprint
548 C Partition local interactions
549 #ifdef MPI
550       call int_bounds(nres-2,loc_start,loc_end)
551       loc_start=loc_start+1
552       loc_end=loc_end+1
553       call int_bounds(nres-2,ithet_start,ithet_end)
554       ithet_start=ithet_start+2
555       ithet_end=ithet_end+2
556       call int_bounds(nct-nnt-2,iturn3_start,iturn3_end) 
557       iturn3_start=iturn3_start+nnt
558       iphi_start=iturn3_start+2
559       iturn3_end=iturn3_end+nnt
560       iphi_end=iturn3_end+2
561       iturn3_start=iturn3_start-1
562       iturn3_end=iturn3_end-1
563       call int_bounds(nres-3,iphi1_start,iphi1_end)
564       iphi1_start=iphi1_start+3
565       iphi1_end=iphi1_end+3
566       call int_bounds(nct-nnt-3,iturn4_start,iturn4_end) 
567       iturn4_start=iturn4_start+nnt
568       iphid_start=iturn4_start+2
569       iturn4_end=iturn4_end+nnt
570       iphid_end=iturn4_end+2
571       iturn4_start=iturn4_start-1
572       iturn4_end=iturn4_end-1
573       call int_bounds(nres-2,ibond_start,ibond_end) 
574       ibond_start=ibond_start+1
575       ibond_end=ibond_end+1
576       call int_bounds(nct-nnt,ibondp_start,ibondp_end) 
577       ibondp_start=ibondp_start+nnt
578       ibondp_end=ibondp_end+nnt
579       call int_bounds1(nres-1,ivec_start,ivec_end) 
580 c      print *,"Processor",myrank,fg_rank,fg_rank1,
581 c     &  " ivec_start",ivec_start," ivec_end",ivec_end
582       iset_start=loc_start+2
583       iset_end=loc_end+2
584       if (ndih_constr.eq.0) then
585         idihconstr_start=1
586         idihconstr_end=0
587       else
588         call int_bounds(ndih_constr,idihconstr_start,idihconstr_end)
589       endif
590 c      nsumgrad=(nres-nnt)*(nres-nnt+1)/2
591 c      nlen=nres-nnt+1
592       nsumgrad=(nres-nnt)*(nres-nnt+1)/2
593       nlen=nres-nnt+1
594       call int_bounds(nsumgrad,ngrad_start,ngrad_end)
595       igrad_start=((2*nlen+1)
596      &    -sqrt(float((2*nlen-1)**2-8*(ngrad_start-1))))/2
597       jgrad_start(igrad_start)=
598      &    ngrad_start-(2*nlen-igrad_start)*(igrad_start-1)/2
599      &    +igrad_start
600       jgrad_end(igrad_start)=nres
601       igrad_end=((2*nlen+1)
602      &    -sqrt(float((2*nlen-1)**2-8*(ngrad_end-1))))/2
603       if (igrad_end.gt.igrad_start) jgrad_start(igrad_end)=igrad_end+1
604       jgrad_end(igrad_end)=ngrad_end-(2*nlen-igrad_end)*(igrad_end-1)/2
605      &    +igrad_end
606       do i=igrad_start+1,igrad_end-1
607         jgrad_start(i)=i+1
608         jgrad_end(i)=nres
609       enddo
610       if (lprint) then 
611         write (*,*) 'Processor:',fg_rank,' CG group',kolor,
612      & ' absolute rank',myrank,
613      & ' loc_start',loc_start,' loc_end',loc_end,
614      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
615      & ' iphi_start',iphi_start,' iphi_end',iphi_end,
616      & ' iphid_start',iphid_start,' iphid_end',iphid_end,
617      & ' ibond_start',ibond_start,' ibond_end',ibond_end,
618      & ' ibondp_start',ibondp_start,' ibondp_end',ibondp_end,
619      & ' iturn3_start',iturn3_start,' iturn3_end',iturn3_end,
620      & ' iturn4_start',iturn4_start,' iturn4_end',iturn4_end,
621      & ' ivec_start',ivec_start,' ivec_end',ivec_end,
622      & ' iset_start',iset_start,' iset_end',iset_end,
623      & ' idihconstr_start',idihconstr_start,' idihconstr_end',
624      &   idihconstr_end
625        write (*,*) 'Processor:',fg_rank,myrank,' igrad_start',
626      &   igrad_start,' igrad_end',igrad_end,' ngrad_start',ngrad_start,
627      &   ' ngrad_end',ngrad_end
628        do i=igrad_start,igrad_end
629          write(*,*) 'Processor:',fg_rank,myrank,i,
630      &    jgrad_start(i),jgrad_end(i)
631        enddo
632       endif
633       if (nfgtasks.gt.1) then
634         call MPI_Allgather(ivec_start,1,MPI_INTEGER,ivec_displ(0),1,
635      &    MPI_INTEGER,FG_COMM1,IERROR)
636         iaux=ivec_end-ivec_start+1
637         call MPI_Allgather(iaux,1,MPI_INTEGER,ivec_count(0),1,
638      &    MPI_INTEGER,FG_COMM1,IERROR)
639         call MPI_Allgather(iset_start-2,1,MPI_INTEGER,iset_displ(0),1,
640      &    MPI_INTEGER,FG_COMM,IERROR)
641         iaux=iset_end-iset_start+1
642         call MPI_Allgather(iaux,1,MPI_INTEGER,iset_count(0),1,
643      &    MPI_INTEGER,FG_COMM,IERROR)
644         call MPI_Allgather(ibond_start,1,MPI_INTEGER,ibond_displ(0),1,
645      &    MPI_INTEGER,FG_COMM,IERROR)
646         iaux=ibond_end-ibond_start+1
647         call MPI_Allgather(iaux,1,MPI_INTEGER,ibond_count(0),1,
648      &    MPI_INTEGER,FG_COMM,IERROR)
649         call MPI_Allgather(ithet_start,1,MPI_INTEGER,ithet_displ(0),1,
650      &    MPI_INTEGER,FG_COMM,IERROR)
651         iaux=ithet_end-ithet_start+1
652         call MPI_Allgather(iaux,1,MPI_INTEGER,ithet_count(0),1,
653      &    MPI_INTEGER,FG_COMM,IERROR)
654         call MPI_Allgather(iphi_start,1,MPI_INTEGER,iphi_displ(0),1,
655      &    MPI_INTEGER,FG_COMM,IERROR)
656         iaux=iphi_end-iphi_start+1
657         call MPI_Allgather(iaux,1,MPI_INTEGER,iphi_count(0),1,
658      &    MPI_INTEGER,FG_COMM,IERROR)
659         call MPI_Allgather(iphi1_start,1,MPI_INTEGER,iphi1_displ(0),1,
660      &    MPI_INTEGER,FG_COMM,IERROR)
661         iaux=iphi1_end-iphi1_start+1
662         call MPI_Allgather(iaux,1,MPI_INTEGER,iphi1_count(0),1,
663      &    MPI_INTEGER,FG_COMM,IERROR)
664         do i=0,maxprocs-1
665           do j=1,maxres
666             ielstart_all(j,i)=0
667             ielend_all(j,i)=0
668           enddo
669         enddo
670         call MPI_Allgather(iturn3_start,1,MPI_INTEGER,
671      &    iturn3_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
672         call MPI_Allgather(iturn4_start,1,MPI_INTEGER,
673      &    iturn4_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
674         call MPI_Allgather(iturn3_end,1,MPI_INTEGER,
675      &    iturn3_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
676         call MPI_Allgather(iturn4_end,1,MPI_INTEGER,
677      &    iturn4_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
678         call MPI_Allgather(iatel_s,1,MPI_INTEGER,
679      &    iatel_s_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
680         call MPI_Allgather(iatel_e,1,MPI_INTEGER,
681      &    iatel_e_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
682         call MPI_Allgather(ielstart(1),maxres,MPI_INTEGER,
683      &    ielstart_all(1,0),maxres,MPI_INTEGER,FG_COMM,IERROR)
684         call MPI_Allgather(ielend(1),maxres,MPI_INTEGER,
685      &    ielend_all(1,0),maxres,MPI_INTEGER,FG_COMM,IERROR)
686         if (lprint) then
687         write (iout,*) "iatel_s_all",(iatel_s_all(i),i=0,nfgtasks)
688         write (iout,*) "iatel_e_all",(iatel_e_all(i),i=0,nfgtasks)
689         write (iout,*) "iturn3_start_all",
690      &    (iturn3_start_all(i),i=0,nfgtasks-1)
691         write (iout,*) "iturn3_end_all",
692      &    (iturn3_end_all(i),i=0,nfgtasks-1)
693         write (iout,*) "iturn4_start_all",
694      &    (iturn4_start_all(i),i=0,nfgtasks-1)
695         write (iout,*) "iturn4_end_all",
696      &    (iturn4_end_all(i),i=0,nfgtasks-1)
697         write (iout,*) "The ielstart_all array"
698         do i=nnt,nct
699           write (iout,'(20i4)') i,(ielstart_all(i,j),j=0,nfgtasks-1)
700         enddo
701         write (iout,*) "The ielend_all array"
702         do i=nnt,nct
703           write (iout,'(20i4)') i,(ielend_all(i,j),j=0,nfgtasks-1)
704         enddo
705         call flush(iout)
706         endif
707         ntask_cont_from=0
708         ntask_cont_to=0
709         itask_cont_from(0)=fg_rank
710         itask_cont_to(0)=fg_rank
711         flag=.false.
712         do ii=iturn3_start,iturn3_end
713           call add_int(ii,ii+2,iturn3_sent(1,ii),
714      &                 ntask_cont_to,itask_cont_to,flag)
715         enddo
716         do ii=iturn4_start,iturn4_end
717           call add_int(ii,ii+3,iturn4_sent(1,ii),
718      &                 ntask_cont_to,itask_cont_to,flag)
719         enddo
720         do ii=iturn3_start,iturn3_end
721           call add_int_from(ii,ii+2,ntask_cont_from,itask_cont_from)
722         enddo
723         do ii=iturn4_start,iturn4_end
724           call add_int_from(ii,ii+3,ntask_cont_from,itask_cont_from)
725         enddo
726         if (lprint) then
727         write (iout,*) "After turn3 ntask_cont_from",ntask_cont_from,
728      &   " ntask_cont_to",ntask_cont_to
729         write (iout,*) "itask_cont_from",
730      &    (itask_cont_from(i),i=1,ntask_cont_from)
731         write (iout,*) "itask_cont_to",
732      &    (itask_cont_to(i),i=1,ntask_cont_to)
733         call flush(iout)
734         endif
735 c        write (iout,*) "Loop forward"
736 c        call flush(iout)
737         do i=iatel_s,iatel_e
738 c          write (iout,*) "from loop i=",i
739 c          call flush(iout)
740           do j=ielstart(i),ielend(i)
741             call add_int_from(i,j,ntask_cont_from,itask_cont_from)
742           enddo
743         enddo
744 c        write (iout,*) "Loop backward iatel_e-1",iatel_e-1,
745 c     &     " iatel_e",iatel_e
746 c        call flush(iout)
747         nat_sent=0
748         do i=iatel_s,iatel_e
749 c          write (iout,*) "i",i," ielstart",ielstart(i),
750 c     &      " ielend",ielend(i)
751 c          call flush(iout)
752           flag=.false.
753           do j=ielstart(i),ielend(i)
754             call add_int(i,j,iint_sent(1,j,nat_sent+1),ntask_cont_to,
755      &                   itask_cont_to,flag)
756           enddo
757           if (flag) then
758             nat_sent=nat_sent+1
759             iat_sent(nat_sent)=i
760           endif
761         enddo
762         if (lprint) then
763         write (iout,*)"After longrange ntask_cont_from",ntask_cont_from,
764      &   " ntask_cont_to",ntask_cont_to
765         write (iout,*) "itask_cont_from",
766      &    (itask_cont_from(i),i=1,ntask_cont_from)
767         write (iout,*) "itask_cont_to",
768      &    (itask_cont_to(i),i=1,ntask_cont_to)
769         call flush(iout)
770         write (iout,*) "iint_sent"
771         do i=1,nat_sent
772           ii=iat_sent(i)
773           write (iout,'(20i4)') ii,(j,(iint_sent(k,j,i),k=1,4),
774      &      j=ielstart(ii),ielend(ii))
775         enddo
776         write (iout,*) "iturn3_sent iturn3_start",iturn3_start,
777      &    " iturn3_end",iturn3_end
778         write (iout,'(20i4)') (i,(iturn3_sent(j,i),j=1,4),
779      &      i=iturn3_start,iturn3_end)
780         write (iout,*) "iturn4_sent iturn4_start",iturn4_start,
781      &    " iturn4_end",iturn4_end
782         write (iout,'(20i4)') (i,(iturn4_sent(j,i),j=1,4),
783      &      i=iturn4_start,iturn4_end)
784         call flush(iout)
785         endif
786         call MPI_Gather(ntask_cont_from,1,MPI_INTEGER,
787      &   ntask_cont_from_all,1,MPI_INTEGER,king,FG_COMM,IERR)
788 c        write (iout,*) "Gather ntask_cont_from ended"
789 c        call flush(iout)
790         call MPI_Gather(itask_cont_from(0),max_fg_procs,MPI_INTEGER,
791      &   itask_cont_from_all(0,0),max_fg_procs,MPI_INTEGER,king,
792      &   FG_COMM,IERR)
793 c        write (iout,*) "Gather itask_cont_from ended"
794 c        call flush(iout)
795         call MPI_Gather(ntask_cont_to,1,MPI_INTEGER,ntask_cont_to_all,
796      &   1,MPI_INTEGER,king,FG_COMM,IERR)
797 c        write (iout,*) "Gather ntask_cont_to ended"
798 c        call flush(iout)
799         call MPI_Gather(itask_cont_to,max_fg_procs,MPI_INTEGER,
800      &   itask_cont_to_all,max_fg_procs,MPI_INTEGER,king,FG_COMM,IERR)
801 c        write (iout,*) "Gather itask_cont_to ended"
802 c        call flush(iout)
803         if (fg_rank.eq.king) then
804           write (iout,*)"Contact receive task map (proc, #tasks, tasks)"
805           do i=0,nfgtasks-1
806             write (iout,'(20i4)') i,ntask_cont_from_all(i),
807      &       (itask_cont_from_all(j,i),j=1,ntask_cont_from_all(i)) 
808           enddo
809           write (iout,*)
810           call flush(iout)
811           write (iout,*) "Contact send task map (proc, #tasks, tasks)"
812           do i=0,nfgtasks-1
813             write (iout,'(20i4)') i,ntask_cont_to_all(i),
814      &       (itask_cont_to_all(j,i),j=1,ntask_cont_to_all(i)) 
815           enddo
816           write (iout,*)
817           call flush(iout)
818 C Check if every send will have a matching receive
819           ncheck_to=0
820           ncheck_from=0
821           do i=0,nfgtasks-1
822             ncheck_to=ncheck_to+ntask_cont_to_all(i)
823             ncheck_from=ncheck_from+ntask_cont_from_all(i)
824           enddo
825           write (iout,*) "Control sums",ncheck_from,ncheck_to
826           if (ncheck_from.ne.ncheck_to) then
827             write (iout,*) "Error: #receive differs from #send."
828             write (iout,*) "Terminating program...!"
829             call flush(iout)
830             flag=.false.
831           else
832             flag=.true.
833             do i=0,nfgtasks-1
834               do j=1,ntask_cont_to_all(i)
835                 ii=itask_cont_to_all(j,i)
836                 do k=1,ntask_cont_from_all(ii)
837                   if (itask_cont_from_all(k,ii).eq.i) then
838                     if(lprint)write(iout,*)"Matching send/receive",i,ii
839                     exit
840                   endif
841                 enddo
842                 if (k.eq.ntask_cont_from_all(ii)+1) then
843                   flag=.false.
844                   write (iout,*) "Error: send by",j," to",ii,
845      &            " would have no matching receive"
846                 endif
847               enddo
848             enddo
849           endif
850           if (.not.flag) then
851             write (iout,*) "Unmatched sends; terminating program"
852             call flush(iout)
853           endif
854         endif
855         call MPI_Bcast(flag,1,MPI_LOGICAL,king,FG_COMM,IERROR)
856 c        write (iout,*) "flag broadcast ended flag=",flag
857 c        call flush(iout)
858         if (.not.flag) then
859           call MPI_Finalize(IERROR)
860           stop "Error in INIT_INT_TABLE: unmatched send/receive."
861         endif
862         call MPI_Comm_group(FG_COMM,fg_group,IERR)
863 c        write (iout,*) "MPI_Comm_group ended"
864 c        call flush(iout)
865         call MPI_Group_incl(fg_group,ntask_cont_from+1,
866      &    itask_cont_from(0),CONT_FROM_GROUP,IERR)
867         call MPI_Group_incl(fg_group,ntask_cont_to+1,itask_cont_to(0),
868      &    CONT_TO_GROUP,IERR)
869         do i=1,nat_sent
870           ii=iat_sent(i)
871           iaux=4*(ielend(ii)-ielstart(ii)+1)
872           call MPI_Group_translate_ranks(fg_group,iaux,
873      &      iint_sent(1,ielstart(ii),i),CONT_TO_GROUP, 
874      &      iint_sent_local(1,ielstart(ii),i),IERR )
875 c          write (iout,*) "Ranks translated i=",i
876 c          call flush(iout)
877         enddo
878         iaux=4*(iturn3_end-iturn3_start+1)
879         call MPI_Group_translate_ranks(fg_group,iaux,
880      &      iturn3_sent(1,iturn3_start),CONT_TO_GROUP,
881      &      iturn3_sent_local(1,iturn3_start),IERR)
882         iaux=4*(iturn4_end-iturn4_start+1)
883         call MPI_Group_translate_ranks(fg_group,iaux,
884      &      iturn4_sent(1,iturn4_start),CONT_TO_GROUP,
885      &      iturn4_sent_local(1,iturn4_start),IERR)
886         if (lprint) then
887         write (iout,*) "iint_sent_local"
888         do i=1,nat_sent
889           ii=iat_sent(i)
890           write (iout,'(20i4)') ii,(j,(iint_sent_local(k,j,i),k=1,4),
891      &      j=ielstart(ii),ielend(ii))
892           call flush(iout)
893         enddo
894         write (iout,*) "iturn3_sent_local iturn3_start",iturn3_start,
895      &    " iturn3_end",iturn3_end
896         write (iout,'(20i4)') (i,(iturn3_sent_local(j,i),j=1,4),
897      &      i=iturn3_start,iturn3_end)
898         write (iout,*) "iturn4_sent_local iturn4_start",iturn4_start,
899      &    " iturn4_end",iturn4_end
900         write (iout,'(20i4)') (i,(iturn4_sent_local(j,i),j=1,4),
901      &      i=iturn4_start,iturn4_end)
902         call flush(iout)
903         endif
904         call MPI_Group_free(fg_group,ierr)
905         call MPI_Group_free(cont_from_group,ierr)
906         call MPI_Group_free(cont_to_group,ierr)
907         call MPI_Type_contiguous(3,MPI_DOUBLE_PRECISION,MPI_UYZ,IERROR)
908         call MPI_Type_commit(MPI_UYZ,IERROR)
909         call MPI_Type_contiguous(18,MPI_DOUBLE_PRECISION,MPI_UYZGRAD,
910      &    IERROR)
911         call MPI_Type_commit(MPI_UYZGRAD,IERROR)
912         call MPI_Type_contiguous(2,MPI_DOUBLE_PRECISION,MPI_MU,IERROR)
913         call MPI_Type_commit(MPI_MU,IERROR)
914         call MPI_Type_contiguous(4,MPI_DOUBLE_PRECISION,MPI_MAT1,IERROR)
915         call MPI_Type_commit(MPI_MAT1,IERROR)
916         call MPI_Type_contiguous(8,MPI_DOUBLE_PRECISION,MPI_MAT2,IERROR)
917         call MPI_Type_commit(MPI_MAT2,IERROR)
918         call MPI_Type_contiguous(6,MPI_DOUBLE_PRECISION,MPI_THET,IERROR)
919         call MPI_Type_commit(MPI_THET,IERROR)
920         call MPI_Type_contiguous(9,MPI_DOUBLE_PRECISION,MPI_GAM,IERROR)
921         call MPI_Type_commit(MPI_GAM,IERROR)
922 #ifndef MATGATHER
923 c 9/22/08 Derived types to send matrices which appear in correlation terms
924         do i=0,nfgtasks-1
925           if (ivec_count(i).eq.ivec_count(0)) then
926             lentyp(i)=0
927           else
928             lentyp(i)=1
929           endif
930         enddo
931         do ind_typ=lentyp(0),lentyp(nfgtasks-1)
932         if (ind_typ.eq.0) then
933           ichunk=ivec_count(0)
934         else
935           ichunk=ivec_count(1)
936         endif
937 c        do i=1,4
938 c          blocklengths(i)=4
939 c        enddo
940 c        displs(1)=0
941 c        do i=2,4
942 c          displs(i)=displs(i-1)+blocklengths(i-1)*maxres
943 c        enddo
944 c        do i=1,4
945 c          blocklengths(i)=blocklengths(i)*ichunk
946 c        enddo
947 c        write (iout,*) "blocklengths and displs"
948 c        do i=1,4
949 c          write (iout,*) i,blocklengths(i),displs(i)
950 c        enddo
951 c        call flush(iout)
952 c        call MPI_Type_indexed(4,blocklengths(1),displs(1),
953 c     &    MPI_DOUBLE_PRECISION,MPI_ROTAT1(ind_typ),IERROR)
954 c        call MPI_Type_commit(MPI_ROTAT1(ind_typ),IERROR)
955 c        write (iout,*) "MPI_ROTAT1",MPI_ROTAT1 
956 c        do i=1,4
957 c          blocklengths(i)=2
958 c        enddo
959 c        displs(1)=0
960 c        do i=2,4
961 c          displs(i)=displs(i-1)+blocklengths(i-1)*maxres
962 c        enddo
963 c        do i=1,4
964 c          blocklengths(i)=blocklengths(i)*ichunk
965 c        enddo
966 c        write (iout,*) "blocklengths and displs"
967 c        do i=1,4
968 c          write (iout,*) i,blocklengths(i),displs(i)
969 c        enddo
970 c        call flush(iout)
971 c        call MPI_Type_indexed(4,blocklengths(1),displs(1),
972 c     &    MPI_DOUBLE_PRECISION,MPI_ROTAT2(ind_typ),IERROR)
973 c        call MPI_Type_commit(MPI_ROTAT2(ind_typ),IERROR)
974 c        write (iout,*) "MPI_ROTAT2",MPI_ROTAT2 
975         do i=1,8
976           blocklengths(i)=2
977         enddo
978         displs(1)=0
979         do i=2,8
980           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
981         enddo
982         do i=1,15
983           blocklengths(i)=blocklengths(i)*ichunk
984         enddo
985         call MPI_Type_indexed(8,blocklengths,displs,
986      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP11(ind_typ),IERROR)
987         call MPI_Type_commit(MPI_PRECOMP11(ind_typ),IERROR)
988         do i=1,8
989           blocklengths(i)=4
990         enddo
991         displs(1)=0
992         do i=2,8
993           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
994         enddo
995         do i=1,15
996           blocklengths(i)=blocklengths(i)*ichunk
997         enddo
998         call MPI_Type_indexed(8,blocklengths,displs,
999      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP12(ind_typ),IERROR)
1000         call MPI_Type_commit(MPI_PRECOMP12(ind_typ),IERROR)
1001         do i=1,6
1002           blocklengths(i)=4
1003         enddo
1004         displs(1)=0
1005         do i=2,6
1006           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1007         enddo
1008         do i=1,6
1009           blocklengths(i)=blocklengths(i)*ichunk
1010         enddo
1011         call MPI_Type_indexed(6,blocklengths,displs,
1012      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP22(ind_typ),IERROR)
1013         call MPI_Type_commit(MPI_PRECOMP22(ind_typ),IERROR)
1014         do i=1,2
1015           blocklengths(i)=8
1016         enddo
1017         displs(1)=0
1018         do i=2,2
1019           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1020         enddo
1021         do i=1,2
1022           blocklengths(i)=blocklengths(i)*ichunk
1023         enddo
1024         call MPI_Type_indexed(2,blocklengths,displs,
1025      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP23(ind_typ),IERROR)
1026         call MPI_Type_commit(MPI_PRECOMP23(ind_typ),IERROR)
1027         do i=1,4
1028           blocklengths(i)=1
1029         enddo
1030         displs(1)=0
1031         do i=2,4
1032           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1033         enddo
1034         do i=1,4
1035           blocklengths(i)=blocklengths(i)*ichunk
1036         enddo
1037         call MPI_Type_indexed(4,blocklengths,displs,
1038      &    MPI_DOUBLE_PRECISION,MPI_ROTAT_OLD(ind_typ),IERROR)
1039         call MPI_Type_commit(MPI_ROTAT_OLD(ind_typ),IERROR)
1040         enddo
1041 #endif
1042       endif
1043       iint_start=ivec_start+1
1044       iint_end=ivec_end+1
1045       do i=0,nfgtasks-1
1046           iint_count(i)=ivec_count(i)
1047           iint_displ(i)=ivec_displ(i)
1048           ivec_displ(i)=ivec_displ(i)-1
1049           iset_displ(i)=iset_displ(i)-1
1050           ithet_displ(i)=ithet_displ(i)-1
1051           iphi_displ(i)=iphi_displ(i)-1
1052           iphi1_displ(i)=iphi1_displ(i)-1
1053           ibond_displ(i)=ibond_displ(i)-1
1054       enddo
1055       if (nfgtasks.gt.1 .and. fg_rank.eq.king 
1056      &     .and. (me.eq.0 .or. .not. out1file)) then
1057         write (iout,*) "IVEC_DISPL, IVEC_COUNT, ISET_START, ISET_COUNT"
1058         do i=0,nfgtasks-1
1059           write (iout,*) i,ivec_displ(i),ivec_count(i),iset_displ(i),
1060      &      iset_count(i)
1061         enddo
1062         write (iout,*) "iphi_start",iphi_start," iphi_end",iphi_end,
1063      &    " iphi1_start",iphi1_start," iphi1_end",iphi1_end
1064         write (iout,*)"IPHI_COUNT, IPHI_DISPL, IPHI1_COUNT, IPHI1_DISPL"
1065         do i=0,nfgtasks-1
1066           write (iout,*) i,iphi_count(i),iphi_displ(i),iphi1_count(i),
1067      &      iphi1_displ(i)
1068         enddo
1069         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
1070      & nele_int_tot,' electrostatic and ',nscp_int_tot,
1071      & ' SC-p interactions','were distributed among',nfgtasks,
1072      & ' fine-grain processors.'
1073       endif
1074 #else
1075       loc_start=2
1076       loc_end=nres-1
1077       ithet_start=3 
1078       ithet_end=nres
1079       iturn3_start=nnt
1080       iturn3_end=nct-3
1081       iturn4_start=nnt
1082       iturn4_end=nct-4
1083       iphi_start=nnt+3
1084       iphi_end=nct
1085       iphi1_start=4
1086       iphi1_end=nres
1087       idihconstr_start=1
1088       idihconstr_end=ndih_constr
1089       iphid_start=iphi_start
1090       iphid_end=iphi_end-1
1091       ibond_start=2
1092       ibond_end=nres-1
1093       ibondp_start=nnt
1094       ibondp_end=nct-1
1095       ivec_start=1
1096       ivec_end=nres-1
1097       iset_start=3
1098       iset_end=nres+1
1099       iint_start=2
1100       iint_end=nres-1
1101 #endif
1102       return
1103       end 
1104 #ifdef MPI
1105 c---------------------------------------------------------------------------
1106       subroutine add_int(ii,jj,itask,ntask_cont_to,itask_cont_to,flag)
1107       implicit none
1108       include "DIMENSIONS"
1109       include "COMMON.INTERACT"
1110       include "COMMON.SETUP"
1111       include "COMMON.IOUNITS"
1112       integer ii,jj,itask(4),ntask_cont_to,
1113      &itask_cont_to(0:max_fg_procs-1)
1114       logical flag
1115       integer iturn3_start_all,iturn3_end_all,iturn4_start_all,
1116      & iturn4_end_all,iatel_s_all,iatel_e_all,ielstart_all,ielend_all
1117       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
1118      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
1119      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
1120      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
1121      & ielend_all(maxres,0:max_fg_procs-1)
1122       integer iproc,isent,k,l
1123 c Determines whether to send interaction ii,jj to other processors; a given
1124 c interaction can be sent to at most 2 processors.
1125 c Sets flag=.true. if interaction ii,jj needs to be sent to at least 
1126 c one processor, otherwise flag is unchanged from the input value.
1127       isent=0
1128       itask(1)=fg_rank
1129       itask(2)=fg_rank
1130       itask(3)=fg_rank
1131       itask(4)=fg_rank
1132 c      write (iout,*) "ii",ii," jj",jj
1133 c Loop over processors to check if anybody could need interaction ii,jj
1134       do iproc=0,fg_rank-1
1135 c Check if the interaction matches any turn3 at iproc
1136         do k=iturn3_start_all(iproc),iturn3_end_all(iproc)
1137           l=k+2
1138           if (k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1
1139      &   .or. k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1 .and. l.eq.jj-1)
1140      &    then 
1141 c            write (iout,*) "turn3 to iproc",iproc," ij",ii,jj,"kl",k,l
1142 c            call flush(iout)
1143             flag=.true.
1144             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1145      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1146               isent=isent+1
1147               itask(isent)=iproc
1148               call add_task(iproc,ntask_cont_to,itask_cont_to)
1149             endif
1150           endif
1151         enddo
1152 C Check if the interaction matches any turn4 at iproc
1153         do k=iturn4_start_all(iproc),iturn4_end_all(iproc)
1154           l=k+3
1155           if (k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1
1156      &   .or. k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1 .and. l.eq.jj-1)
1157      &    then 
1158 c            write (iout,*) "turn3 to iproc",iproc," ij",ii,jj," kl",k,l
1159 c            call flush(iout)
1160             flag=.true.
1161             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1162      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1163               isent=isent+1
1164               itask(isent)=iproc
1165               call add_task(iproc,ntask_cont_to,itask_cont_to)
1166             endif
1167           endif
1168         enddo
1169         if (iatel_s_all(iproc).gt.0 .and. iatel_e_all(iproc).gt.0 .and. 
1170      &  iatel_s_all(iproc).le.ii-1 .and. iatel_e_all(iproc).ge.ii-1)then
1171           if (ielstart_all(ii-1,iproc).le.jj-1.and.
1172      &        ielend_all(ii-1,iproc).ge.jj-1) then
1173             flag=.true.
1174             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1175      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1176               isent=isent+1
1177               itask(isent)=iproc
1178               call add_task(iproc,ntask_cont_to,itask_cont_to)
1179             endif
1180           endif
1181           if (ielstart_all(ii-1,iproc).le.jj+1.and.
1182      &        ielend_all(ii-1,iproc).ge.jj+1) then
1183             flag=.true.
1184             if (iproc.ne.itask(1).and.iproc.ne.itask(2)
1185      &        .and.iproc.ne.itask(3).and.iproc.ne.itask(4)) then
1186               isent=isent+1
1187               itask(isent)=iproc
1188               call add_task(iproc,ntask_cont_to,itask_cont_to)
1189             endif
1190           endif
1191         endif
1192       enddo
1193       return
1194       end
1195 c---------------------------------------------------------------------------
1196       subroutine add_int_from(ii,jj,ntask_cont_from,itask_cont_from)
1197       implicit none
1198       include "DIMENSIONS"
1199       include "COMMON.INTERACT"
1200       include "COMMON.SETUP"
1201       include "COMMON.IOUNITS"
1202       integer ii,jj,itask(2),ntask_cont_from,
1203      & itask_cont_from(0:max_fg_procs-1)
1204       logical flag
1205       integer iturn3_start_all,iturn3_end_all,iturn4_start_all,
1206      & iturn4_end_all,iatel_s_all,iatel_e_all,ielstart_all,ielend_all
1207       common /przechowalnia/ iturn3_start_all(0:max_fg_procs),
1208      & iturn3_end_all(0:max_fg_procs),iturn4_start_all(0:max_fg_procs),
1209      & iturn4_end_all(0:max_fg_procs),iatel_s_all(0:max_fg_procs),
1210      &iatel_e_all(0:max_fg_procs),ielstart_all(maxres,0:max_fg_procs-1),
1211      & ielend_all(maxres,0:max_fg_procs-1)
1212       integer iproc,k,l
1213       do iproc=fg_rank+1,nfgtasks-1
1214         do k=iturn3_start_all(iproc),iturn3_end_all(iproc)
1215           l=k+2
1216           if (k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1.and.l.eq.jj-1 
1217      &   .or. k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1) 
1218      &    then
1219 c            write (iout,*)"turn3 from iproc",iproc," ij",ii,jj," kl",k,l
1220             call add_task(iproc,ntask_cont_from,itask_cont_from)
1221           endif
1222         enddo 
1223         do k=iturn4_start_all(iproc),iturn4_end_all(iproc)
1224           l=k+3
1225           if (k.eq.ii+1 .and. l.eq.jj+1 .or. k.eq.ii+1.and.l.eq.jj-1 
1226      &   .or. k.eq.ii-1 .and. l.eq.jj-1 .or. k.eq.ii-1 .and. l.eq.jj+1) 
1227      &    then
1228 c            write (iout,*)"turn4 from iproc",iproc," ij",ii,jj," kl",k,l
1229             call add_task(iproc,ntask_cont_from,itask_cont_from)
1230           endif
1231         enddo 
1232         if (iatel_s_all(iproc).gt.0 .and. iatel_e_all(iproc).gt.0) then
1233           if (ii+1.ge.iatel_s_all(iproc).and.ii+1.le.iatel_e_all(iproc))
1234      &    then
1235             if (jj+1.ge.ielstart_all(ii+1,iproc).and.
1236      &          jj+1.le.ielend_all(ii+1,iproc)) then
1237               call add_task(iproc,ntask_cont_from,itask_cont_from)
1238             endif            
1239             if (jj-1.ge.ielstart_all(ii+1,iproc).and.
1240      &          jj-1.le.ielend_all(ii+1,iproc)) then
1241               call add_task(iproc,ntask_cont_from,itask_cont_from)
1242             endif
1243           endif
1244           if (ii-1.ge.iatel_s_all(iproc).and.ii-1.le.iatel_e_all(iproc))
1245      &    then
1246             if (jj-1.ge.ielstart_all(ii-1,iproc).and.
1247      &          jj-1.le.ielend_all(ii-1,iproc)) then
1248               call add_task(iproc,ntask_cont_from,itask_cont_from)
1249             endif
1250             if (jj+1.ge.ielstart_all(ii-1,iproc).and.
1251      &          jj+1.le.ielend_all(ii-1,iproc)) then
1252                call add_task(iproc,ntask_cont_from,itask_cont_from)
1253             endif
1254           endif
1255         endif
1256       enddo
1257       return
1258       end
1259 c---------------------------------------------------------------------------
1260       subroutine add_task(iproc,ntask_cont,itask_cont)
1261       implicit none
1262       include "DIMENSIONS"
1263       integer iproc,ntask_cont,itask_cont(0:max_fg_procs-1)
1264       integer ii
1265       do ii=1,ntask_cont
1266         if (itask_cont(ii).eq.iproc) return
1267       enddo
1268       ntask_cont=ntask_cont+1
1269       itask_cont(ntask_cont)=iproc
1270       return
1271       end
1272 c---------------------------------------------------------------------------
1273       subroutine int_bounds(total_ints,lower_bound,upper_bound)
1274       implicit real*8 (a-h,o-z)
1275       include 'DIMENSIONS'
1276       include 'mpif.h'
1277       include 'COMMON.SETUP'
1278       integer total_ints,lower_bound,upper_bound
1279       integer int4proc(0:max_fg_procs),sint4proc(0:max_fg_procs)
1280       nint=total_ints/nfgtasks
1281       do i=1,nfgtasks
1282         int4proc(i-1)=nint
1283       enddo
1284       nexcess=total_ints-nint*nfgtasks
1285       do i=1,nexcess
1286         int4proc(nfgtasks-i)=int4proc(nfgtasks-i)+1
1287       enddo
1288       lower_bound=0
1289       do i=0,fg_rank-1
1290         lower_bound=lower_bound+int4proc(i)
1291       enddo 
1292       upper_bound=lower_bound+int4proc(fg_rank)
1293       lower_bound=lower_bound+1
1294       return
1295       end
1296 c---------------------------------------------------------------------------
1297       subroutine int_bounds1(total_ints,lower_bound,upper_bound)
1298       implicit real*8 (a-h,o-z)
1299       include 'DIMENSIONS'
1300       include 'mpif.h'
1301       include 'COMMON.SETUP'
1302       integer total_ints,lower_bound,upper_bound
1303       integer int4proc(0:max_fg_procs),sint4proc(0:max_fg_procs)
1304       nint=total_ints/nfgtasks1
1305       do i=1,nfgtasks1
1306         int4proc(i-1)=nint
1307       enddo
1308       nexcess=total_ints-nint*nfgtasks1
1309       do i=1,nexcess
1310         int4proc(nfgtasks1-i)=int4proc(nfgtasks1-i)+1
1311       enddo
1312       lower_bound=0
1313       do i=0,fg_rank1-1
1314         lower_bound=lower_bound+int4proc(i)
1315       enddo 
1316       upper_bound=lower_bound+int4proc(fg_rank1)
1317       lower_bound=lower_bound+1
1318       return
1319       end
1320 c---------------------------------------------------------------------------
1321       subroutine int_partition(int_index,lower_index,upper_index,atom,
1322      & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
1323       implicit real*8 (a-h,o-z)
1324       include 'DIMENSIONS'
1325       include 'COMMON.IOUNITS'
1326       integer int_index,lower_index,upper_index,atom,at_start,at_end,
1327      & first_atom,last_atom,int_gr,jat_start,jat_end
1328       logical lprn
1329       lprn=.false.
1330       if (lprn) write (iout,*) 'int_index=',int_index
1331       int_index_old=int_index
1332       int_index=int_index+last_atom-first_atom+1
1333       if (lprn) 
1334      &   write (iout,*) 'int_index=',int_index,
1335      &               ' int_index_old',int_index_old,
1336      &               ' lower_index=',lower_index,
1337      &               ' upper_index=',upper_index,
1338      &               ' atom=',atom,' first_atom=',first_atom,
1339      &               ' last_atom=',last_atom
1340       if (int_index.ge.lower_index) then
1341         int_gr=int_gr+1
1342         if (at_start.eq.0) then
1343           at_start=atom
1344           jat_start=first_atom-1+lower_index-int_index_old
1345         else
1346           jat_start=first_atom
1347         endif
1348         if (lprn) write (iout,*) 'jat_start',jat_start
1349         if (int_index.ge.upper_index) then
1350           at_end=atom
1351           jat_end=first_atom-1+upper_index-int_index_old
1352           return1
1353         else
1354           jat_end=last_atom
1355         endif
1356         if (lprn) write (iout,*) 'jat_end',jat_end
1357       endif
1358       return
1359       end
1360 #endif
1361 c------------------------------------------------------------------------------
1362       subroutine hpb_partition
1363       implicit real*8 (a-h,o-z)
1364       include 'DIMENSIONS'
1365 #ifdef MPI
1366       include 'mpif.h'
1367 #endif
1368       include 'COMMON.SBRIDGE'
1369       include 'COMMON.IOUNITS'
1370       include 'COMMON.SETUP'
1371 #ifdef MPI
1372       call int_bounds(nhpb,link_start,link_end)
1373       write (iout,*) 'Processor',fg_rank,' CG group',kolor,
1374      &  ' absolute rank',MyRank,
1375      &  ' nhpb',nhpb,' link_start=',link_start,
1376      &  ' link_end',link_end
1377 #else
1378       link_start=1
1379       link_end=nhpb
1380 #endif
1381       return
1382       end