source/cluster/wham/src-M/initialize_p.F

   1       subroutine initialize
   2 C
   3 C Define constants and zero out tables.
   4 C
   5       implicit real*8 (a-h,o-z)
   6       include 'DIMENSIONS'
   7       include 'sizesclu.dat'
   8 #ifdef MPI
   9       include 'mpif.h'
  10 #endif
  11       include 'COMMON.IOUNITS'
  12       include 'COMMON.CHAIN'
  13       include 'COMMON.INTERACT'
  14       include 'COMMON.GEO'
  15       include 'COMMON.LOCAL'
  16       include 'COMMON.TORSION'
  17       include 'COMMON.FFIELD'
  18       include 'COMMON.SBRIDGE'
  19       include 'COMMON.MINIM'
  20       include 'COMMON.DERIV'
  21       include "COMMON.NAMES"
  22       include "COMMON.TIME1"
  23 C
  24 C The following is just to define auxiliary variables used in angle conversion
  25 C
  26       pi=4.0D0*datan(1.0D0)
  27       dwapi=2.0D0*pi
  28       dwapi3=dwapi/3.0D0
  29       pipol=0.5D0*pi
  30       deg2rad=pi/180.0D0
  31       rad2deg=1.0D0/deg2rad
  32       angmin=10.0D0*deg2rad
  33       Rgas = 1.987D-3
  34 C
  35 C Define I/O units.
  36 C
  37       inp=    1
  38       iout=   2
  39       ipdbin= 3
  40       ipdb=   7
  41       imol2= 18
  42       jplot= 19
  43       jstatin=10
  44       imol2=  4
  45       igeom=  8
  46       intin=  9
  47       ithep= 11
  48       irotam=12
  49       itorp= 13
  50       itordp= 23
  51       ielep= 14
  52       isidep=15
  53       isidep1=22
  54       iscpp=25
  55       icbase=16
  56       ifourier=20
  57       istat= 17
  58       ibond=28
  59       isccor=29
  60       jrms=30
  61 C
  62 C Set default weights of the energy terms.
  63 C
  64       wlong=1.0D0
  65       welec=1.0D0
  66       wtor =1.0D0
  67       wang =1.0D0
  68       wscloc=1.0D0
  69       wstrain=1.0D0
  70 C
  71 C Zero out tables.
  72 C
  73       ndih_constr=0
  74       do i=1,maxres2
  75         do j=1,3
  76           c(j,i)=0.0D0
  77           dc(j,i)=0.0D0
  78         enddo
  79       enddo
  80       do i=1,maxres
  81         do j=1,3
  82           xloc(j,i)=0.0D0
  83         enddo
  84       enddo
  85       do i=1,ntyp
  86         do j=1,ntyp
  87           aa(i,j)=0.0D0
  88           bb(i,j)=0.0D0
  89           augm(i,j)=0.0D0
  90           sigma(i,j)=0.0D0
  91           r0(i,j)=0.0D0
  92           chi(i,j)=0.0D0
  93         enddo
  94         do j=1,2
  95           bad(i,j)=0.0D0
  96         enddo
  97         chip(i)=0.0D0
  98         alp(i)=0.0D0
  99         sigma0(i)=0.0D0
 100         sigii(i)=0.0D0
 101         rr0(i)=0.0D0
 102         a0thet(i)=0.0D0
 103         do j=1,2
 104           athet(j,i)=0.0D0
 105           bthet(j,i)=0.0D0
 106         enddo
 107         do j=0,3
 108           polthet(j,i)=0.0D0
 109         enddo
 110         do j=1,3
 111           gthet(j,i)=0.0D0
 112         enddo
 113         theta0(i)=0.0D0
 114         sig0(i)=0.0D0
 115         sigc0(i)=0.0D0
 116         do j=1,maxlob
 117           bsc(j,i)=0.0D0
 118           do k=1,3
 119             censc(k,j,i)=0.0D0
 120           enddo
 121           do k=1,3
 122             do l=1,3
 123               gaussc(l,k,j,i)=0.0D0
 124             enddo
 125           enddo
 126           nlob(i)=0
 127         enddo
 128       enddo
 129       nlob(ntyp1)=0
 130       dsc(ntyp1)=0.0D0
 131       do i=1,maxtor
 132         itortyp(i)=0
 133         do j=1,maxtor
 134           do k=1,maxterm
 135             v1(k,j,i)=0.0D0
 136             v2(k,j,i)=0.0D0
 137           enddo
 138         enddo
 139       enddo
 140       do i=1,maxres
 141         itype(i)=0
 142         itel(i)=0
 143       enddo
 144 C Initialize the bridge arrays
 145       ns=0
 146       nss=0
 147       nhpb=0
 148       do i=1,maxss
 149         iss(i)=0
 150       enddo
 151       do i=1,maxss
 152         dhpb(i)=0.0D0
 153       enddo
 154       do i=1,maxss
 155         ihpb(i)=0
 156         jhpb(i)=0
 157       enddo
 158 C
 159 C Initialize timing.
 160 C
 161       call set_timers
 162 C
 163 C Initialize variables used in minimization.
 164 C
 165 c     maxfun=5000
 166 c     maxit=2000
 167       maxfun=500
 168       maxit=200
 169       tolf=1.0D-2
 170       rtolf=5.0D-4
 171 C
 172 C Initialize the variables responsible for the mode of gradient storage.
 173 C
 174       nfl=0
 175       icg=1
 176       do i=1,14
 177         do j=1,14
 178           if (print_order(i).eq.j) then
 179             iw(print_order(i))=j
 180             goto 1121
 181           endif
 182         enddo
 183 1121    continue
 184       enddo
 185       calc_grad=.false.
 186 C Set timers and counters for the respective routines
 187       t_func = 0.0d0
 188       t_grad = 0.0d0
 189       t_fhel = 0.0d0
 190       t_fbet = 0.0d0
 191       t_ghel = 0.0d0
 192       t_gbet = 0.0d0
 193       t_viol = 0.0d0
 194       t_gviol = 0.0d0
 195       n_func = 0
 196       n_grad = 0
 197       n_fhel = 0
 198       n_fbet = 0
 199       n_ghel = 0
 200       n_gbet = 0
 201       n_viol = 0
 202       n_gviol = 0
 203       n_map = 0
 204 #ifndef SPLITELE
 205       nprint_ene=nprint_ene-1
 206 #endif
 207       return
 208       end
 209 c-------------------------------------------------------------------------
 210       block data nazwy
 211       implicit real*8 (a-h,o-z)
 212       include 'DIMENSIONS'
 213       include 'sizesclu.dat'
 214       include 'COMMON.NAMES'
 215       include 'COMMON.FFIELD'
 216       data restyp /
 217      &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
 218      &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','D'/
 219       data onelet /
 220      &'C','M','F','I','L','V','W','Y','A','G','T',
 221      &'S','Q','N','E','D','H','R','K','P','X'/
 222       data potname /'LJ','LJK','BP','GB','GBV'/
 223       data ename /
 224      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
 225      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
 226      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP",
 227      &   "ESTR","ESCCOR","EVDW2_14","","EVDW_T"/
 228       data wname /
 229      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
 230      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
 231      &   "WHPB","WVDWPP","WBOND","WSCCOR","WSCP14","","WSC"/
 232       data nprint_ene /21/
 233       data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,
 234      &  16,15,17,20,21/
 235       end
 236 c---------------------------------------------------------------------------
 237       subroutine init_int_table
 238       implicit real*8 (a-h,o-z)
 239       include 'DIMENSIONS'
 240       include 'sizesclu.dat'
 241 #ifdef MPI
 242       include 'mpif.h'
 243 #endif
 244 #ifdef MPL
 245       include 'COMMON.INFO'
 246 #endif
 247       include 'COMMON.CHAIN'
 248       include 'COMMON.INTERACT'
 249       include 'COMMON.LOCAL'
 250       include 'COMMON.SBRIDGE'
 251       include 'COMMON.IOUNITS'
 252       logical scheck,lprint
 253 #ifdef MPL
 254       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
 255      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
 256 C... Determine the numbers of start and end SC-SC interaction
 257 C... to deal with by current processor.
 258       lprint=.false.
 259       if (lprint)
 260      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
 261       n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
 262       MyRank=MyID-(MyGroup-1)*fgProcs
 263       call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
 264       if (lprint)
 265      &  write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 266      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
 267      &  ' my_sc_inde',my_sc_inde
 268       ind_sctint=0
 269       iatsc_s=0
 270       iatsc_e=0
 271 #endif
 272       lprint=.false.
 273       do i=1,maxres
 274         nint_gr(i)=0
 275         nscp_gr(i)=0
 276         do j=1,maxint_gr
 277           istart(i,1)=0
 278           iend(i,1)=0
 279           ielstart(i)=0
 280           ielend(i)=0
 281           iscpstart(i,1)=0
 282           iscpend(i,1)=0
 283         enddo
 284       enddo
 285       ind_scint=0
 286       ind_scint_old=0
 287 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
 288 cd   &   (ihpb(i),jhpb(i),i=1,nss)
 289       do i=nnt,nct-1
 290         scheck=.false.
 291         do ii=1,nss
 292           if (ihpb(ii).eq.i+nres) then
 293             scheck=.true.
 294             jj=jhpb(ii)-nres
 295             goto 10
 296           endif
 297         enddo
 298    10   continue
 299 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
 300         if (scheck) then
 301           if (jj.eq.i+1) then
 302 #ifdef MPL
 303             write (iout,*) 'jj=i+1'
 304             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 305      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 306 #else
 307             nint_gr(i)=1
 308             istart(i,1)=i+2
 309             iend(i,1)=nct
 310 #endif
 311           else if (jj.eq.nct) then
 312 #ifdef MPL
 313             write (iout,*) 'jj=nct'
 314             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 315      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 316 #else
 317             nint_gr(i)=1
 318             istart(i,1)=i+1
 319             iend(i,1)=nct-1
 320 #endif
 321           else
 322 #ifdef MPL
 323             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 324      & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 325             ii=nint_gr(i)+1
 326             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 327      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
 328 #else
 329             nint_gr(i)=2
 330             istart(i,1)=i+1
 331             iend(i,1)=jj-1
 332             istart(i,2)=jj+1
 333             iend(i,2)=nct
 334 #endif
 335           endif
 336         else
 337 #ifdef MPL
 338           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 339      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 340 #else
 341           nint_gr(i)=1
 342           istart(i,1)=i+1
 343           iend(i,1)=nct
 344           ind_scint=int_scint+nct-i
 345 #endif
 346         endif
 347 #ifdef MPL
 348         ind_scint_old=ind_scint
 349 #endif
 350       enddo
 351    12 continue
 352 #ifndef MPL
 353       iatsc_s=nnt
 354       iatsc_e=nct-1
 355 #endif
 356 #ifdef MPL
 357       if (lprint) then
 358         write (iout,*) 'Processor',MyID,' Group',MyGroup
 359         write (iout,*) 'iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
 360       endif
 361 #endif
 362       if (lprint) then
 363       write (iout,'(a)') 'Interaction array:'
 364       do i=iatsc_s,iatsc_e
 365         write (iout,'(i3,2(2x,2i3))')
 366      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
 367       enddo
 368       endif
 369       ispp=2
 370 #ifdef MPL
 371 C Now partition the electrostatic-interaction array
 372       npept=nct-nnt
 373       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
 374       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
 375       if (lprint)
 376      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 377      &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
 378      &               ' my_ele_inde',my_ele_inde
 379       iatel_s=0
 380       iatel_e=0
 381       ind_eleint=0
 382       ind_eleint_old=0
 383       do i=nnt,nct-3
 384         ijunk=0
 385         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
 386      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
 387       enddo ! i
 388    13 continue
 389 #else
 390       iatel_s=nnt
 391       iatel_e=nct-3
 392       do i=iatel_s,iatel_e
 393         ielstart(i)=i+2
 394         ielend(i)=nct-1
 395       enddo
 396 #endif
 397       if (lprint) then
 398         write (iout,'(a)') 'Electrostatic interaction array:'
 399         do i=iatel_s,iatel_e
 400           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
 401         enddo
 402       endif ! lprint
 403 c     iscp=3
 404       iscp=2
 405 C Partition the SC-p interaction array
 406 #ifdef MPL
 407       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
 408       call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
 409       if (lprint)
 410      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 411      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
 412      &               ' my_scp_inde',my_scp_inde
 413       iatscp_s=0
 414       iatscp_e=0
 415       ind_scpint=0
 416       ind_scpint_old=0
 417       do i=nnt,nct-1
 418         if (i.lt.nnt+iscp) then
 419 cd        write (iout,*) 'i.le.nnt+iscp'
 420           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 421      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
 422      &      iscpend(i,1),*14)
 423         else if (i.gt.nct-iscp) then
 424 cd        write (iout,*) 'i.gt.nct-iscp'
 425           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 426      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 427      &      iscpend(i,1),*14)
 428         else
 429           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 430      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 431      &      iscpend(i,1),*14)
 432           ii=nscp_gr(i)+1
 433           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 434      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
 435      &      iscpend(i,ii),*14)
 436         endif
 437       enddo ! i
 438    14 continue
 439 #else
 440       iatscp_s=nnt
 441       iatscp_e=nct-1
 442       do i=nnt,nct-1
 443         if (i.lt.nnt+iscp) then
 444           nscp_gr(i)=1
 445           iscpstart(i,1)=i+iscp
 446           iscpend(i,1)=nct
 447         elseif (i.gt.nct-iscp) then
 448           nscp_gr(i)=1
 449           iscpstart(i,1)=nnt
 450           iscpend(i,1)=i-iscp
 451         else
 452           nscp_gr(i)=2
 453           iscpstart(i,1)=nnt
 454           iscpend(i,1)=i-iscp
 455           iscpstart(i,2)=i+iscp
 456           iscpend(i,2)=nct
 457         endif
 458       enddo ! i
 459 #endif
 460       if (lprint) then
 461         write (iout,'(a)') 'SC-p interaction array:'
 462         do i=iatscp_s,iatscp_e
 463           write (iout,'(i3,2(2x,2i3))')
 464      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
 465         enddo
 466       endif ! lprint
 467 C Partition local interactions
 468 #ifdef MPL
 469       call int_bounds(nres-2,loc_start,loc_end)
 470       loc_start=loc_start+1
 471       loc_end=loc_end+1
 472       call int_bounds(nres-2,ithet_start,ithet_end)
 473       ithet_start=ithet_start+2
 474       ithet_end=ithet_end+2
 475       call int_bounds(nct-nnt-2,iphi_start,iphi_end)
 476       iphi_start=iphi_start+nnt+2
 477       iphi_end=iphi_end+nnt+2
 478       if (lprint) then
 479         write (iout,*) 'Processor:',MyID,
 480      & ' loc_start',loc_start,' loc_end',loc_end,
 481      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 482      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 483         write (*,*) 'Processor:',MyID,
 484      & ' loc_start',loc_start,' loc_end',loc_end,
 485      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 486      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 487       endif
 488       if (fgprocs.gt.1 .and. MyID.eq.BossID) then
 489         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
 490      & nele_int_tot,' electrostatic and ',nscp_int_tot,
 491      & ' SC-p interactions','were distributed among',fgprocs,
 492      & ' fine-grain processors.'
 493       endif
 494 #else
 495       loc_start=2
 496       loc_end=nres-1
 497       ithet_start=3
 498       ithet_end=nres
 499       iphi_start=nnt+3
 500       iphi_end=nct
 501 #endif
 502       return
 503       end
 504 c---------------------------------------------------------------------------
 505       subroutine int_partition(int_index,lower_index,upper_index,atom,
 506      & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
 507       implicit real*8 (a-h,o-z)
 508       include 'DIMENSIONS'
 509       include 'COMMON.IOUNITS'
 510       integer int_index,lower_index,upper_index,atom,at_start,at_end,
 511      & first_atom,last_atom,int_gr,jat_start,jat_end
 512       logical lprn
 513       lprn=.false.
 514       if (lprn) write (iout,*) 'int_index=',int_index
 515       int_index_old=int_index
 516       int_index=int_index+last_atom-first_atom+1
 517       if (lprn)
 518      &   write (iout,*) 'int_index=',int_index,
 519      &               ' int_index_old',int_index_old,
 520      &               ' lower_index=',lower_index,
 521      &               ' upper_index=',upper_index,
 522      &               ' atom=',atom,' first_atom=',first_atom,
 523      &               ' last_atom=',last_atom
 524       if (int_index.ge.lower_index) then
 525         int_gr=int_gr+1
 526         if (at_start.eq.0) then
 527           at_start=atom
 528           jat_start=first_atom-1+lower_index-int_index_old
 529         else
 530           jat_start=first_atom
 531         endif
 532         if (lprn) write (iout,*) 'jat_start',jat_start
 533         if (int_index.ge.upper_index) then
 534           at_end=atom
 535           jat_end=first_atom-1+upper_index-int_index_old
 536           return1
 537         else
 538           jat_end=last_atom
 539         endif
 540         if (lprn) write (iout,*) 'jat_end',jat_end
 541       endif
 542       return
 543       end
 544 c------------------------------------------------------------------------------
 545       subroutine hpb_partition
 546       implicit real*8 (a-h,o-z)
 547       include 'DIMENSIONS'
 548       include 'COMMON.SBRIDGE'
 549       include 'COMMON.IOUNITS'
 550 #ifdef MPL
 551       include 'COMMON.INFO'
 552       call int_bounds(nhpb,link_start,link_end)
 553 #else
 554       link_start=1
 555       link_end=nhpb
 556 #endif
 557 cd    write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 558 cd   &  ' nhpb',nhpb,' link_start=',link_start,
 559 cd   &  ' link_end',link_end
 560       return
 561       end