62706fba9a6bf6b405ea9f0e6fb5cb10c173d441
[unres.git] / source / unres / src_MD / 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       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           athet(j,i)=0.0D0
168           bthet(j,i)=0.0D0
169         enddo
170         do j=0,3
171           polthet(j,i)=0.0D0
172         enddo
173         do j=1,3
174           gthet(j,i)=0.0D0
175         enddo
176         theta0(i)=0.0D0
177         sig0(i)=0.0D0
178         sigc0(i)=0.0D0
179         do j=1,maxlob
180           bsc(j,i)=0.0D0
181           do k=1,3
182             censc(k,j,i)=0.0D0
183           enddo
184           do k=1,3
185             do l=1,3
186               gaussc(l,k,j,i)=0.0D0
187             enddo
188           enddo
189           nlob(i)=0
190         enddo
191       enddo
192       nlob(ntyp1)=0
193       dsc(ntyp1)=0.0D0
194       do i=-maxtor,maxtor
195         itortyp(i)=0
196        do iblock=1,2
197         do j=-maxtor,maxtor
198           do k=1,maxterm
199             v1(k,j,i,iblock)=0.0D0
200             v2(k,j,i,iblock)=0.0D0
201           enddo
202         enddo
203         enddo
204       enddo
205       do iblock=1,2
206        do i=-maxtor,maxtor
207         do j=-maxtor,maxtor
208          do k=-maxtor,maxtor
209           do l=1,maxtermd_1
210             v1c(1,l,i,j,k,iblock)=0.0D0
211             v1s(1,l,i,j,k,iblock)=0.0D0
212             v1c(2,l,i,j,k,iblock)=0.0D0
213             v1s(2,l,i,j,k,iblock)=0.0D0
214           enddo !l
215           do l=1,maxtermd_2
216            do m=1,maxtermd_2
217             v2c(m,l,i,j,k,iblock)=0.0D0
218             v2s(m,l,i,j,k,iblock)=0.0D0
219            enddo !m
220           enddo !l
221         enddo !k
222        enddo !j
223       enddo !i
224       enddo !iblock
225       do i=1,maxres
226         itype(i)=0
227         itel(i)=0
228       enddo
229 C Initialize the bridge arrays
230       ns=0
231       nss=0 
232       nhpb=0
233       do i=1,maxss
234         iss(i)=0
235       enddo
236       do i=1,maxdim
237         dhpb(i)=0.0D0
238       enddo
239       do i=1,maxres
240         ihpb(i)=0
241         jhpb(i)=0
242       enddo
243 C
244 C Initialize timing.
245 C
246       call set_timers
247 C
248 C Initialize variables used in minimization.
249 C   
250 c     maxfun=5000
251 c     maxit=2000
252       maxfun=500
253       maxit=200
254       tolf=1.0D-2
255       rtolf=5.0D-4
256
257 C Initialize the variables responsible for the mode of gradient storage.
258 C
259       nfl=0
260       icg=1
261 C
262 C Initialize constants used to split the energy into long- and short-range
263 C components
264 C
265       r_cut=2.0d0
266       rlamb=0.3d0
267 #ifndef SPLITELE
268       nprint_ene=nprint_ene-1
269 #endif
270       return
271       end
272 c-------------------------------------------------------------------------
273       block data nazwy
274       implicit real*8 (a-h,o-z)
275       include 'DIMENSIONS'
276       include 'COMMON.NAMES'
277       include 'COMMON.FFIELD'
278       data restyp /
279      &'DD' ,'DPR','DLY','DAR','DHI','DAS','DGL','DSG','DGN','DSN','DTH',
280      &'DYY','DAL','DTY','DTR','DVA','DLE','DIL','DPN','MED','DCY','ZER',
281      &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
282      &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','D'/
283       data onelet /
284      &'z','p','k','r','h','d','e','n','q','s','t','g',
285      &'a','y','w','v','l','i','f','m','c','x',
286      &'C','M','F','I','L','V','W','Y','A','G','T',
287      &'S','Q','N','E','D','H','R','K','P','X'/
288       data potname /'LJ','LJK','BP','GB','GBV'/
289       data ename /
290      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
291      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
292      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB ","EVDWPP ",
293      &   "ESTR ","EVDW2_14 ","UCONST ", "      ","ESCCOR"," "," "/
294       data wname /
295      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
296      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
297      &   "WSTRAIN","WVDWPP","WBOND","SCAL14","     ","    ","WSCCOR",
298      &   " "," "/
299       data nprint_ene /20/
300       data print_order/1,2,3,11,12,13,14,4,5,6,7,8,9,10,19,18,15,17,16,
301      & 21,0,0,0/
302       end 
303 c---------------------------------------------------------------------------
304       subroutine init_int_table
305       implicit real*8 (a-h,o-z)
306       include 'DIMENSIONS'
307 #ifdef MPI
308       include 'mpif.h'
309       integer blocklengths(15),displs(15)
310 #endif
311       include 'COMMON.CONTROL'
312       include 'COMMON.SETUP'
313       include 'COMMON.CHAIN'
314       include 'COMMON.INTERACT'
315       include 'COMMON.LOCAL'
316       include 'COMMON.SBRIDGE'
317       include 'COMMON.TORCNSTR'
318       include 'COMMON.IOUNITS'
319       include 'COMMON.DERIV'
320       include 'COMMON.CONTACTS'
321       common /przechowalnia/ iturn3_start_all(0:MaxProcs),
322      & iturn3_end_all(0:MaxProcs),iturn4_start_all(0:MaxProcs),
323      & iturn4_end_all(0:MaxProcs),iatel_s_all(0:MaxProcs),
324      & iatel_e_all(0:MaxProcs),ielstart_all(maxres,0:MaxProcs-1),
325      & ielend_all(maxres,0:MaxProcs-1),
326      & ntask_cont_from_all(0:max_fg_procs-1),
327      & itask_cont_from_all(0:max_fg_procs-1,0:max_fg_procs-1),
328      & ntask_cont_to_all(0:max_fg_procs-1),
329      & itask_cont_to_all(0:max_fg_procs-1,0:max_fg_procs-1)
330       integer FG_GROUP,CONT_FROM_GROUP,CONT_TO_GROUP
331       logical scheck,lprint,flag
332 #ifdef MPI
333       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
334      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
335 C... Determine the numbers of start and end SC-SC interaction 
336 C... to deal with by current processor.
337       do i=0,nfgtasks-1
338         itask_cont_from(i)=fg_rank
339         itask_cont_to(i)=fg_rank
340       enddo
341       lprint=.false.
342       if (lprint)
343      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
344       n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
345       call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
346       if (lprint)
347      &  write (iout,*) 'Processor',fg_rank,' CG group',kolor,
348      &  ' absolute rank',MyRank,
349      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
350      &  ' my_sc_inde',my_sc_inde
351       ind_sctint=0
352       iatsc_s=0
353       iatsc_e=0
354 #endif
355 c      lprint=.false.
356       do i=1,maxres
357         nint_gr(i)=0
358         nscp_gr(i)=0
359         do j=1,maxint_gr
360           istart(i,1)=0
361           iend(i,1)=0
362           ielstart(i)=0
363           ielend(i)=0
364           iscpstart(i,1)=0
365           iscpend(i,1)=0    
366         enddo
367       enddo
368       ind_scint=0
369       ind_scint_old=0
370 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
371 cd   &   (ihpb(i),jhpb(i),i=1,nss)
372       do i=nnt,nct-1
373         scheck=.false.
374         do ii=1,nss
375           if (ihpb(ii).eq.i+nres) then
376             scheck=.true.
377             jj=jhpb(ii)-nres
378             goto 10
379           endif
380         enddo
381    10   continue
382 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
383         if (scheck) then
384           if (jj.eq.i+1) then
385 #ifdef MPI
386 c            write (iout,*) 'jj=i+1'
387             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
388      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
389 #else
390             nint_gr(i)=1
391             istart(i,1)=i+2
392             iend(i,1)=nct
393 #endif
394           else if (jj.eq.nct) then
395 #ifdef MPI
396 c            write (iout,*) 'jj=nct'
397             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
398      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
399 #else
400             nint_gr(i)=1
401             istart(i,1)=i+1
402             iend(i,1)=nct-1
403 #endif
404           else
405 #ifdef MPI
406             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
407      & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
408             ii=nint_gr(i)+1
409             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
410      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
411 #else
412             nint_gr(i)=2
413             istart(i,1)=i+1
414             iend(i,1)=jj-1
415             istart(i,2)=jj+1
416             iend(i,2)=nct
417 #endif
418           endif
419         else
420 #ifdef MPI
421           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
422      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
423 #else
424           nint_gr(i)=1
425           istart(i,1)=i+1
426           iend(i,1)=nct
427           ind_scint=ind_scint+nct-i
428 #endif
429         endif
430 #ifdef MPI
431         ind_scint_old=ind_scint
432 #endif
433       enddo
434    12 continue
435 #ifndef MPI
436       iatsc_s=nnt
437       iatsc_e=nct-1
438 #endif
439 #ifdef MPI
440       if (lprint) write (*,*) 'Processor',fg_rank,' CG Group',kolor,
441      &   ' absolute rank',myrank,' iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
442 #endif
443       if (lprint) then
444       write (iout,'(a)') 'Interaction array:'
445       do i=iatsc_s,iatsc_e
446         write (iout,'(i3,2(2x,2i3))') 
447      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
448       enddo
449       endif
450       ispp=4
451 #ifdef MPI
452 C Now partition the electrostatic-interaction array
453       npept=nct-nnt
454       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
455       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
456       if (lprint)
457      & write (*,*) 'Processor',fg_rank,' CG group',kolor,
458      &  ' absolute rank',MyRank,
459      &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
460      &               ' my_ele_inde',my_ele_inde
461       iatel_s=0
462       iatel_e=0
463       ind_eleint=0
464       ind_eleint_old=0
465       do i=nnt,nct-3
466         ijunk=0
467         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
468      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
469       enddo ! i 
470    13 continue
471       if (iatel_s.eq.0) iatel_s=1
472       nele_int_tot_vdw=(npept-2)*(npept-2+1)/2
473 c      write (iout,*) "nele_int_tot_vdw",nele_int_tot_vdw
474       call int_bounds(nele_int_tot_vdw,my_ele_inds_vdw,my_ele_inde_vdw)
475 c      write (iout,*) "my_ele_inds_vdw",my_ele_inds_vdw,
476 c     & " my_ele_inde_vdw",my_ele_inde_vdw
477       ind_eleint_vdw=0
478       ind_eleint_vdw_old=0
479       iatel_s_vdw=0
480       iatel_e_vdw=0
481       do i=nnt,nct-3
482         ijunk=0
483         call int_partition(ind_eleint_vdw,my_ele_inds_vdw,
484      &    my_ele_inde_vdw,i,
485      &    iatel_s_vdw,iatel_e_vdw,i+2,nct-1,ijunk,ielstart_vdw(i),
486      &    ielend_vdw(i),*15)
487 c        write (iout,*) i," ielstart_vdw",ielstart_vdw(i),
488 c     &   " ielend_vdw",ielend_vdw(i)
489       enddo ! i 
490       if (iatel_s_vdw.eq.0) iatel_s_vdw=1
491    15 continue
492 #else
493       iatel_s=nnt
494       iatel_e=nct-5
495       do i=iatel_s,iatel_e
496         ielstart(i)=i+4
497         ielend(i)=nct-1
498       enddo
499       iatel_s_vdw=nnt
500       iatel_e_vdw=nct-3
501       do i=iatel_s_vdw,iatel_e_vdw
502         ielstart_vdw(i)=i+2
503         ielend_vdw(i)=nct-1
504       enddo
505 #endif
506       if (lprint) then
507         write (*,'(a)') 'Processor',fg_rank,' CG group',kolor,
508      &  ' absolute rank',MyRank
509         write (iout,*) 'Electrostatic interaction array:'
510         do i=iatel_s,iatel_e
511           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
512         enddo
513       endif ! lprint
514 c     iscp=3
515       iscp=2
516 C Partition the SC-p interaction array
517 #ifdef MPI
518       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
519       call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
520       if (lprint) write (iout,*) 'Processor',fg_rank,' CG group',kolor,
521      &  ' absolute rank',myrank,
522      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
523      &               ' my_scp_inde',my_scp_inde
524       iatscp_s=0
525       iatscp_e=0
526       ind_scpint=0
527       ind_scpint_old=0
528       do i=nnt,nct-1
529         if (i.lt.nnt+iscp) then
530 cd        write (iout,*) 'i.le.nnt+iscp'
531           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
532      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
533      &      iscpend(i,1),*14)
534         else if (i.gt.nct-iscp) then
535 cd        write (iout,*) 'i.gt.nct-iscp'
536           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
537      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
538      &      iscpend(i,1),*14)
539         else
540           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
541      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
542      &      iscpend(i,1),*14)
543           ii=nscp_gr(i)+1
544           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
545      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
546      &      iscpend(i,ii),*14)
547         endif
548       enddo ! i
549    14 continue
550 #else
551       iatscp_s=nnt
552       iatscp_e=nct-1
553       do i=nnt,nct-1
554         if (i.lt.nnt+iscp) then
555           nscp_gr(i)=1
556           iscpstart(i,1)=i+iscp
557           iscpend(i,1)=nct
558         elseif (i.gt.nct-iscp) then
559           nscp_gr(i)=1
560           iscpstart(i,1)=nnt
561           iscpend(i,1)=i-iscp
562         else
563           nscp_gr(i)=2
564           iscpstart(i,1)=nnt
565           iscpend(i,1)=i-iscp
566           iscpstart(i,2)=i+iscp
567           iscpend(i,2)=nct
568         endif 
569       enddo ! i
570 #endif
571       if (lprint) then
572         write (iout,'(a)') 'SC-p interaction array:'
573         do i=iatscp_s,iatscp_e
574           write (iout,'(i3,2(2x,2i3))') 
575      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
576         enddo
577       endif ! lprint
578 C Partition local interactions
579 #ifdef MPI
580       call int_bounds(nres-2,loc_start,loc_end)
581       loc_start=loc_start+1
582       loc_end=loc_end+1
583       call int_bounds(nres-2,ithet_start,ithet_end)
584       ithet_start=ithet_start+2
585       ithet_end=ithet_end+2
586       call int_bounds(nct-nnt-2,iturn3_start,iturn3_end) 
587       iturn3_start=iturn3_start+nnt
588       iphi_start=iturn3_start+2
589       iturn3_end=iturn3_end+nnt
590       iphi_end=iturn3_end+2
591       iturn3_start=iturn3_start-1
592       iturn3_end=iturn3_end-1
593       call int_bounds(nres-3,itau_start,itau_end) 
594       itau_start=itau_start+3
595       itau_end=itau_end+3
596       call int_bounds(nres-3,iphi1_start,iphi1_end)
597       iphi1_start=iphi1_start+3
598       iphi1_end=iphi1_end+3
599       call int_bounds(nct-nnt-3,iturn4_start,iturn4_end) 
600       iturn4_start=iturn4_start+nnt
601       iphid_start=iturn4_start+2
602       iturn4_end=iturn4_end+nnt
603       iphid_end=iturn4_end+2
604       iturn4_start=iturn4_start-1
605       iturn4_end=iturn4_end-1
606       call int_bounds(nres-2,ibond_start,ibond_end) 
607       ibond_start=ibond_start+1
608       ibond_end=ibond_end+1
609       call int_bounds(nct-nnt,ibondp_start,ibondp_end) 
610       ibondp_start=ibondp_start+nnt
611       ibondp_end=ibondp_end+nnt
612       call int_bounds1(nres-1,ivec_start,ivec_end) 
613       print *,"Processor",myrank,fg_rank,fg_rank1,
614      &  " ivec_start",ivec_start," ivec_end",ivec_end
615       iset_start=loc_start+2
616       iset_end=loc_end+2
617       if (ndih_constr.eq.0) then
618         idihconstr_start=1
619         idihconstr_end=0
620       else
621         call int_bounds(ndih_constr,idihconstr_start,idihconstr_end)
622       endif
623       nsumgrad=(nres-nnt)*(nres-nnt+1)/2
624       nlen=nres-nnt+1
625       call int_bounds(nsumgrad,ngrad_start,ngrad_end)
626       igrad_start=((2*nlen+1)
627      &    -sqrt(float((2*nlen-1)**2-8*(ngrad_start-1))))/2
628       jgrad_start(igrad_start)=
629      &    ngrad_start-(2*nlen-igrad_start)*(igrad_start-1)/2
630      &    +igrad_start
631       jgrad_end(igrad_start)=nres
632       igrad_end=((2*nlen+1)
633      &    -sqrt(float((2*nlen-1)**2-8*(ngrad_end-1))))/2
634       if (igrad_end.gt.igrad_start) jgrad_start(igrad_end)=igrad_end+1
635       jgrad_end(igrad_end)=ngrad_end-(2*nlen-igrad_end)*(igrad_end-1)/2
636      &    +igrad_end
637       do i=igrad_start+1,igrad_end-1
638         jgrad_start(i)=i+1
639         jgrad_end(i)=nres
640       enddo
641       if (lprint) then 
642         write (*,*) 'Processor:',fg_rank,' CG group',kolor,
643      & ' absolute rank',myrank,
644      & ' loc_start',loc_start,' loc_end',loc_end,
645      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
646      & ' iphi_start',iphi_start,' iphi_end',iphi_end,
647      & ' iphid_start',iphid_start,' iphid_end',iphid_end,
648      & ' ibond_start',ibond_start,' ibond_end',ibond_end,
649      & ' ibondp_start',ibondp_start,' ibondp_end',ibondp_end,
650      & ' iturn3_start',iturn3_start,' iturn3_end',iturn3_end,
651      & ' iturn4_start',iturn4_start,' iturn4_end',iturn4_end,
652      & ' ivec_start',ivec_start,' ivec_end',ivec_end,
653      & ' iset_start',iset_start,' iset_end',iset_end,
654      & ' idihconstr_start',idihconstr_start,' idihconstr_end',
655      &   idihconstr_end
656        write (*,*) 'Processor:',fg_rank,myrank,' igrad_start',
657      &   igrad_start,' igrad_end',igrad_end,' ngrad_start',ngrad_start,
658      &   ' ngrad_end',ngrad_end
659        do i=igrad_start,igrad_end
660          write(*,*) 'Processor:',fg_rank,myrank,i,
661      &    jgrad_start(i),jgrad_end(i)
662        enddo
663       endif
664       if (nfgtasks.gt.1) then
665         call MPI_Allgather(ivec_start,1,MPI_INTEGER,ivec_displ(0),1,
666      &    MPI_INTEGER,FG_COMM1,IERROR)
667         iaux=ivec_end-ivec_start+1
668         call MPI_Allgather(iaux,1,MPI_INTEGER,ivec_count(0),1,
669      &    MPI_INTEGER,FG_COMM1,IERROR)
670         call MPI_Allgather(iset_start-2,1,MPI_INTEGER,iset_displ(0),1,
671      &    MPI_INTEGER,FG_COMM,IERROR)
672         iaux=iset_end-iset_start+1
673         call MPI_Allgather(iaux,1,MPI_INTEGER,iset_count(0),1,
674      &    MPI_INTEGER,FG_COMM,IERROR)
675         call MPI_Allgather(ibond_start,1,MPI_INTEGER,ibond_displ(0),1,
676      &    MPI_INTEGER,FG_COMM,IERROR)
677         iaux=ibond_end-ibond_start+1
678         call MPI_Allgather(iaux,1,MPI_INTEGER,ibond_count(0),1,
679      &    MPI_INTEGER,FG_COMM,IERROR)
680         call MPI_Allgather(ithet_start,1,MPI_INTEGER,ithet_displ(0),1,
681      &    MPI_INTEGER,FG_COMM,IERROR)
682         iaux=ithet_end-ithet_start+1
683         call MPI_Allgather(iaux,1,MPI_INTEGER,ithet_count(0),1,
684      &    MPI_INTEGER,FG_COMM,IERROR)
685         call MPI_Allgather(iphi_start,1,MPI_INTEGER,iphi_displ(0),1,
686      &    MPI_INTEGER,FG_COMM,IERROR)
687         iaux=iphi_end-iphi_start+1
688         call MPI_Allgather(iaux,1,MPI_INTEGER,iphi_count(0),1,
689      &    MPI_INTEGER,FG_COMM,IERROR)
690         call MPI_Allgather(iphi1_start,1,MPI_INTEGER,iphi1_displ(0),1,
691      &    MPI_INTEGER,FG_COMM,IERROR)
692         iaux=iphi1_end-iphi1_start+1
693         call MPI_Allgather(iaux,1,MPI_INTEGER,iphi1_count(0),1,
694      &    MPI_INTEGER,FG_COMM,IERROR)
695         do i=0,maxprocs-1
696           do j=1,maxres
697             ielstart_all(j,i)=0
698             ielend_all(j,i)=0
699           enddo
700         enddo
701         call MPI_Allgather(iturn3_start,1,MPI_INTEGER,
702      &    iturn3_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
703         call MPI_Allgather(iturn4_start,1,MPI_INTEGER,
704      &    iturn4_start_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
705         call MPI_Allgather(iturn3_end,1,MPI_INTEGER,
706      &    iturn3_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
707         call MPI_Allgather(iturn4_end,1,MPI_INTEGER,
708      &    iturn4_end_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
709         call MPI_Allgather(iatel_s,1,MPI_INTEGER,
710      &    iatel_s_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
711         call MPI_Allgather(iatel_e,1,MPI_INTEGER,
712      &    iatel_e_all(0),1,MPI_INTEGER,FG_COMM,IERROR)
713         call MPI_Allgather(ielstart(1),maxres,MPI_INTEGER,
714      &    ielstart_all(1,0),maxres,MPI_INTEGER,FG_COMM,IERROR)
715         call MPI_Allgather(ielend(1),maxres,MPI_INTEGER,
716      &    ielend_all(1,0),maxres,MPI_INTEGER,FG_COMM,IERROR)
717         if (lprint) then
718         write (iout,*) "iatel_s_all",(iatel_s_all(i),i=0,nfgtasks)
719         write (iout,*) "iatel_e_all",(iatel_e_all(i),i=0,nfgtasks)
720         write (iout,*) "iturn3_start_all",
721      &    (iturn3_start_all(i),i=0,nfgtasks-1)
722         write (iout,*) "iturn3_end_all",
723      &    (iturn3_end_all(i),i=0,nfgtasks-1)
724         write (iout,*) "iturn4_start_all",
725      &    (iturn4_start_all(i),i=0,nfgtasks-1)
726         write (iout,*) "iturn4_end_all",
727      &    (iturn4_end_all(i),i=0,nfgtasks-1)
728         write (iout,*) "The ielstart_all array"
729         do i=nnt,nct
730           write (iout,'(20i4)') i,(ielstart_all(i,j),j=0,nfgtasks-1)
731         enddo
732         write (iout,*) "The ielend_all array"
733         do i=nnt,nct
734           write (iout,'(20i4)') i,(ielend_all(i,j),j=0,nfgtasks-1)
735         enddo
736         call flush(iout)
737         endif
738         ntask_cont_from=0
739         ntask_cont_to=0
740         itask_cont_from(0)=fg_rank
741         itask_cont_to(0)=fg_rank
742         flag=.false.
743         do ii=iturn3_start,iturn3_end
744           call add_int(ii,ii+2,iturn3_sent(1,ii),
745      &                 ntask_cont_to,itask_cont_to,flag)
746         enddo
747         do ii=iturn4_start,iturn4_end
748           call add_int(ii,ii+3,iturn4_sent(1,ii),
749      &                 ntask_cont_to,itask_cont_to,flag)
750         enddo
751         do ii=iturn3_start,iturn3_end
752           call add_int_from(ii,ii+2,ntask_cont_from,itask_cont_from)
753         enddo
754         do ii=iturn4_start,iturn4_end
755           call add_int_from(ii,ii+3,ntask_cont_from,itask_cont_from)
756         enddo
757         if (lprint) then
758         write (iout,*) "After turn3 ntask_cont_from",ntask_cont_from,
759      &   " ntask_cont_to",ntask_cont_to
760         write (iout,*) "itask_cont_from",
761      &    (itask_cont_from(i),i=1,ntask_cont_from)
762         write (iout,*) "itask_cont_to",
763      &    (itask_cont_to(i),i=1,ntask_cont_to)
764         call flush(iout)
765         endif
766 c        write (iout,*) "Loop forward"
767 c        call flush(iout)
768         do i=iatel_s,iatel_e
769 c          write (iout,*) "from loop i=",i
770 c          call flush(iout)
771           do j=ielstart(i),ielend(i)
772             call add_int_from(i,j,ntask_cont_from,itask_cont_from)
773           enddo
774         enddo
775 c        write (iout,*) "Loop backward iatel_e-1",iatel_e-1,
776 c     &     " iatel_e",iatel_e
777 c        call flush(iout)
778         nat_sent=0
779         do i=iatel_s,iatel_e
780 c          write (iout,*) "i",i," ielstart",ielstart(i),
781 c     &      " ielend",ielend(i)
782 c          call flush(iout)
783           flag=.false.
784           do j=ielstart(i),ielend(i)
785             call add_int(i,j,iint_sent(1,j,nat_sent+1),ntask_cont_to,
786      &                   itask_cont_to,flag)
787           enddo
788           if (flag) then
789             nat_sent=nat_sent+1
790             iat_sent(nat_sent)=i
791           endif
792         enddo
793         if (lprint) then
794         write (iout,*)"After longrange ntask_cont_from",ntask_cont_from,
795      &   " ntask_cont_to",ntask_cont_to
796         write (iout,*) "itask_cont_from",
797      &    (itask_cont_from(i),i=1,ntask_cont_from)
798         write (iout,*) "itask_cont_to",
799      &    (itask_cont_to(i),i=1,ntask_cont_to)
800         call flush(iout)
801         write (iout,*) "iint_sent"
802         do i=1,nat_sent
803           ii=iat_sent(i)
804           write (iout,'(20i4)') ii,(j,(iint_sent(k,j,i),k=1,4),
805      &      j=ielstart(ii),ielend(ii))
806         enddo
807         write (iout,*) "iturn3_sent iturn3_start",iturn3_start,
808      &    " iturn3_end",iturn3_end
809         write (iout,'(20i4)') (i,(iturn3_sent(j,i),j=1,4),
810      &      i=iturn3_start,iturn3_end)
811         write (iout,*) "iturn4_sent iturn4_start",iturn4_start,
812      &    " iturn4_end",iturn4_end
813         write (iout,'(20i4)') (i,(iturn4_sent(j,i),j=1,4),
814      &      i=iturn4_start,iturn4_end)
815         call flush(iout)
816         endif
817         call MPI_Gather(ntask_cont_from,1,MPI_INTEGER,
818      &   ntask_cont_from_all,1,MPI_INTEGER,king,FG_COMM,IERR)
819 c        write (iout,*) "Gather ntask_cont_from ended"
820 c        call flush(iout)
821         call MPI_Gather(itask_cont_from(0),max_fg_procs,MPI_INTEGER,
822      &   itask_cont_from_all(0,0),max_fg_procs,MPI_INTEGER,king,
823      &   FG_COMM,IERR)
824 c        write (iout,*) "Gather itask_cont_from ended"
825 c        call flush(iout)
826         call MPI_Gather(ntask_cont_to,1,MPI_INTEGER,ntask_cont_to_all,
827      &   1,MPI_INTEGER,king,FG_COMM,IERR)
828 c        write (iout,*) "Gather ntask_cont_to ended"
829 c        call flush(iout)
830         call MPI_Gather(itask_cont_to,max_fg_procs,MPI_INTEGER,
831      &   itask_cont_to_all,max_fg_procs,MPI_INTEGER,king,FG_COMM,IERR)
832 c        write (iout,*) "Gather itask_cont_to ended"
833 c        call flush(iout)
834         if (fg_rank.eq.king) then
835           write (iout,*)"Contact receive task map (proc, #tasks, tasks)"
836           do i=0,nfgtasks-1
837             write (iout,'(20i4)') i,ntask_cont_from_all(i),
838      &       (itask_cont_from_all(j,i),j=1,ntask_cont_from_all(i)) 
839           enddo
840           write (iout,*)
841           call flush(iout)
842           write (iout,*) "Contact send task map (proc, #tasks, tasks)"
843           do i=0,nfgtasks-1
844             write (iout,'(20i4)') i,ntask_cont_to_all(i),
845      &       (itask_cont_to_all(j,i),j=1,ntask_cont_to_all(i)) 
846           enddo
847           write (iout,*)
848           call flush(iout)
849 C Check if every send will have a matching receive
850           ncheck_to=0
851           ncheck_from=0
852           do i=0,nfgtasks-1
853             ncheck_to=ncheck_to+ntask_cont_to_all(i)
854             ncheck_from=ncheck_from+ntask_cont_from_all(i)
855           enddo
856           write (iout,*) "Control sums",ncheck_from,ncheck_to
857           if (ncheck_from.ne.ncheck_to) then
858             write (iout,*) "Error: #receive differs from #send."
859             write (iout,*) "Terminating program...!"
860             call flush(iout)
861             flag=.false.
862           else
863             flag=.true.
864             do i=0,nfgtasks-1
865               do j=1,ntask_cont_to_all(i)
866                 ii=itask_cont_to_all(j,i)
867                 do k=1,ntask_cont_from_all(ii)
868                   if (itask_cont_from_all(k,ii).eq.i) then
869                     if(lprint)write(iout,*)"Matching send/receive",i,ii
870                     exit
871                   endif
872                 enddo
873                 if (k.eq.ntask_cont_from_all(ii)+1) then
874                   flag=.false.
875                   write (iout,*) "Error: send by",j," to",ii,
876      &            " would have no matching receive"
877                 endif
878               enddo
879             enddo
880           endif
881           if (.not.flag) then
882             write (iout,*) "Unmatched sends; terminating program"
883             call flush(iout)
884           endif
885         endif
886         call MPI_Bcast(flag,1,MPI_LOGICAL,king,FG_COMM,IERROR)
887 c        write (iout,*) "flag broadcast ended flag=",flag
888 c        call flush(iout)
889         if (.not.flag) then
890           call MPI_Finalize(IERROR)
891           stop "Error in INIT_INT_TABLE: unmatched send/receive."
892         endif
893         call MPI_Comm_group(FG_COMM,fg_group,IERR)
894 c        write (iout,*) "MPI_Comm_group ended"
895 c        call flush(iout)
896         call MPI_Group_incl(fg_group,ntask_cont_from+1,
897      &    itask_cont_from(0),CONT_FROM_GROUP,IERR)
898         call MPI_Group_incl(fg_group,ntask_cont_to+1,itask_cont_to(0),
899      &    CONT_TO_GROUP,IERR)
900         do i=1,nat_sent
901           ii=iat_sent(i)
902           iaux=4*(ielend(ii)-ielstart(ii)+1)
903           call MPI_Group_translate_ranks(fg_group,iaux,
904      &      iint_sent(1,ielstart(ii),i),CONT_TO_GROUP, 
905      &      iint_sent_local(1,ielstart(ii),i),IERR )
906 c          write (iout,*) "Ranks translated i=",i
907 c          call flush(iout)
908         enddo
909         iaux=4*(iturn3_end-iturn3_start+1)
910         call MPI_Group_translate_ranks(fg_group,iaux,
911      &      iturn3_sent(1,iturn3_start),CONT_TO_GROUP,
912      &      iturn3_sent_local(1,iturn3_start),IERR)
913         iaux=4*(iturn4_end-iturn4_start+1)
914         call MPI_Group_translate_ranks(fg_group,iaux,
915      &      iturn4_sent(1,iturn4_start),CONT_TO_GROUP,
916      &      iturn4_sent_local(1,iturn4_start),IERR)
917         if (lprint) then
918         write (iout,*) "iint_sent_local"
919         do i=1,nat_sent
920           ii=iat_sent(i)
921           write (iout,'(20i4)') ii,(j,(iint_sent_local(k,j,i),k=1,4),
922      &      j=ielstart(ii),ielend(ii))
923           call flush(iout)
924         enddo
925         write (iout,*) "iturn3_sent_local iturn3_start",iturn3_start,
926      &    " iturn3_end",iturn3_end
927         write (iout,'(20i4)') (i,(iturn3_sent_local(j,i),j=1,4),
928      &      i=iturn3_start,iturn3_end)
929         write (iout,*) "iturn4_sent_local iturn4_start",iturn4_start,
930      &    " iturn4_end",iturn4_end
931         write (iout,'(20i4)') (i,(iturn4_sent_local(j,i),j=1,4),
932      &      i=iturn4_start,iturn4_end)
933         call flush(iout)
934         endif
935         call MPI_Group_free(fg_group,ierr)
936         call MPI_Group_free(cont_from_group,ierr)
937         call MPI_Group_free(cont_to_group,ierr)
938         call MPI_Type_contiguous(3,MPI_DOUBLE_PRECISION,MPI_UYZ,IERROR)
939         call MPI_Type_commit(MPI_UYZ,IERROR)
940         call MPI_Type_contiguous(18,MPI_DOUBLE_PRECISION,MPI_UYZGRAD,
941      &    IERROR)
942         call MPI_Type_commit(MPI_UYZGRAD,IERROR)
943         call MPI_Type_contiguous(2,MPI_DOUBLE_PRECISION,MPI_MU,IERROR)
944         call MPI_Type_commit(MPI_MU,IERROR)
945         call MPI_Type_contiguous(4,MPI_DOUBLE_PRECISION,MPI_MAT1,IERROR)
946         call MPI_Type_commit(MPI_MAT1,IERROR)
947         call MPI_Type_contiguous(8,MPI_DOUBLE_PRECISION,MPI_MAT2,IERROR)
948         call MPI_Type_commit(MPI_MAT2,IERROR)
949         call MPI_Type_contiguous(6,MPI_DOUBLE_PRECISION,MPI_THET,IERROR)
950         call MPI_Type_commit(MPI_THET,IERROR)
951         call MPI_Type_contiguous(9,MPI_DOUBLE_PRECISION,MPI_GAM,IERROR)
952         call MPI_Type_commit(MPI_GAM,IERROR)
953 #ifndef MATGATHER
954 c 9/22/08 Derived types to send matrices which appear in correlation terms
955         do i=0,nfgtasks-1
956           if (ivec_count(i).eq.ivec_count(0)) then
957             lentyp(i)=0
958           else
959             lentyp(i)=1
960           endif
961         enddo
962         do ind_typ=lentyp(0),lentyp(nfgtasks-1)
963         if (ind_typ.eq.0) then
964           ichunk=ivec_count(0)
965         else
966           ichunk=ivec_count(1)
967         endif
968 c        do i=1,4
969 c          blocklengths(i)=4
970 c        enddo
971 c        displs(1)=0
972 c        do i=2,4
973 c          displs(i)=displs(i-1)+blocklengths(i-1)*maxres
974 c        enddo
975 c        do i=1,4
976 c          blocklengths(i)=blocklengths(i)*ichunk
977 c        enddo
978 c        write (iout,*) "blocklengths and displs"
979 c        do i=1,4
980 c          write (iout,*) i,blocklengths(i),displs(i)
981 c        enddo
982 c        call flush(iout)
983 c        call MPI_Type_indexed(4,blocklengths(1),displs(1),
984 c     &    MPI_DOUBLE_PRECISION,MPI_ROTAT1(ind_typ),IERROR)
985 c        call MPI_Type_commit(MPI_ROTAT1(ind_typ),IERROR)
986 c        write (iout,*) "MPI_ROTAT1",MPI_ROTAT1 
987 c        do i=1,4
988 c          blocklengths(i)=2
989 c        enddo
990 c        displs(1)=0
991 c        do i=2,4
992 c          displs(i)=displs(i-1)+blocklengths(i-1)*maxres
993 c        enddo
994 c        do i=1,4
995 c          blocklengths(i)=blocklengths(i)*ichunk
996 c        enddo
997 c        write (iout,*) "blocklengths and displs"
998 c        do i=1,4
999 c          write (iout,*) i,blocklengths(i),displs(i)
1000 c        enddo
1001 c        call flush(iout)
1002 c        call MPI_Type_indexed(4,blocklengths(1),displs(1),
1003 c     &    MPI_DOUBLE_PRECISION,MPI_ROTAT2(ind_typ),IERROR)
1004 c        call MPI_Type_commit(MPI_ROTAT2(ind_typ),IERROR)
1005 c        write (iout,*) "MPI_ROTAT2",MPI_ROTAT2 
1006         do i=1,8
1007           blocklengths(i)=2
1008         enddo
1009         displs(1)=0
1010         do i=2,8
1011           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1012         enddo
1013         do i=1,15
1014           blocklengths(i)=blocklengths(i)*ichunk
1015         enddo
1016         call MPI_Type_indexed(8,blocklengths,displs,
1017      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP11(ind_typ),IERROR)
1018         call MPI_Type_commit(MPI_PRECOMP11(ind_typ),IERROR)
1019         do i=1,8
1020           blocklengths(i)=4
1021         enddo
1022         displs(1)=0
1023         do i=2,8
1024           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1025         enddo
1026         do i=1,15
1027           blocklengths(i)=blocklengths(i)*ichunk
1028         enddo
1029         call MPI_Type_indexed(8,blocklengths,displs,
1030      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP12(ind_typ),IERROR)
1031         call MPI_Type_commit(MPI_PRECOMP12(ind_typ),IERROR)
1032         do i=1,6
1033           blocklengths(i)=4
1034         enddo
1035         displs(1)=0
1036         do i=2,6
1037           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1038         enddo
1039         do i=1,6
1040           blocklengths(i)=blocklengths(i)*ichunk
1041         enddo
1042         call MPI_Type_indexed(6,blocklengths,displs,
1043      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP22(ind_typ),IERROR)
1044         call MPI_Type_commit(MPI_PRECOMP22(ind_typ),IERROR)
1045         do i=1,2
1046           blocklengths(i)=8
1047         enddo
1048         displs(1)=0
1049         do i=2,2
1050           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1051         enddo
1052         do i=1,2
1053           blocklengths(i)=blocklengths(i)*ichunk
1054         enddo
1055         call MPI_Type_indexed(2,blocklengths,displs,
1056      &    MPI_DOUBLE_PRECISION,MPI_PRECOMP23(ind_typ),IERROR)
1057         call MPI_Type_commit(MPI_PRECOMP23(ind_typ),IERROR)
1058         do i=1,4
1059           blocklengths(i)=1
1060         enddo
1061         displs(1)=0
1062         do i=2,4
1063           displs(i)=displs(i-1)+blocklengths(i-1)*maxres
1064         enddo
1065         do i=1,4
1066           blocklengths(i)=blocklengths(i)*ichunk
1067         enddo
1068         call MPI_Type_indexed(4,blocklengths,displs,
1069      &    MPI_DOUBLE_PRECISION,MPI_ROTAT_OLD(ind_typ),IERROR)
1070         call MPI_Type_commit(MPI_ROTAT_OLD(ind_typ),IERROR)
1071         enddo
1072 #endif
1073       endif
1074       iint_start=ivec_start+1
1075       iint_end=ivec_end+1
1076       do i=0,nfgtasks-1
1077           iint_count(i)=ivec_count(i)
1078           iint_displ(i)=ivec_displ(i)
1079           ivec_displ(i)=ivec_displ(i)-1
1080           iset_displ(i)=iset_displ(i)-1
1081           ithet_displ(i)=ithet_displ(i)-1
1082           iphi_displ(i)=iphi_displ(i)-1
1083           iphi1_displ(i)=iphi1_displ(i)-1
1084           ibond_displ(i)=ibond_displ(i)-1
1085       enddo
1086       if (nfgtasks.gt.1 .and. fg_rank.eq.king 
1087      &     .and. (me.eq.0 .or. out1file)) then
1088         write (iout,*) "IVEC_DISPL, IVEC_COUNT, ISET_START, ISET_COUNT"
1089         do i=0,nfgtasks-1
1090           write (iout,*) i,ivec_displ(i),ivec_count(i),iset_displ(i),
1091      &      iset_count(i)
1092         enddo
1093         write (iout,*) "iphi_start",iphi_start," iphi_end",iphi_end,
1094      &    " iphi1_start",iphi1_start," iphi1_end",iphi1_end
1095         write (iout,*)"IPHI_COUNT, IPHI_DISPL, IPHI1_COUNT, IPHI1_DISPL"
1096         do i=0,nfgtasks-1
1097           write (iout,*) i,iphi_count(i),iphi_displ(i),iphi1_count(i),
1098      &      iphi1_displ(i)
1099         enddo
1100         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
1101      & nele_int_tot,' electrostatic and ',nscp_int_tot,
1102      & ' SC-p interactions','were distributed among',nfgtasks,
1103      & ' fine-grain processors.'
1104       endif
1105 #else
1106       loc_start=2
1107       loc_end=nres-1
1108       ithet_start=3 
1109       ithet_end=nres
1110       iturn3_start=nnt
1111       iturn3_end=nct-3
1112       iturn4_start=nnt
1113       iturn4_end=nct-4
1114       iphi_start=nnt+3
1115       iphi_end=nct
1116       iphi1_start=4
1117       iphi1_end=nres
1118       idihconstr_start=1
1119       idihconstr_end=ndih_constr
1120       iphid_start=iphi_start
1121       iphid_end=iphi_end-1
1122       itau_start=4
1123       itau_end=nres
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