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