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