source/unres/src_MD-M/initialize_p.F

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