source/cluster/wham/src-NEWSC/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       include 'COMMON.INTERACT'
 217       data restyp /
 218      &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
 219      &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','D'/
 220       data onelet /
 221      &'C','M','F','I','L','V','W','Y','A','G','T',
 222      &'S','Q','N','E','D','H','R','K','P','X'/
 223       data potname /'LJ','LJK','BP','GB','GBV','MM'/
 224       data ename /
 225      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
 226      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
 227      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP",
 228      &   "ESTR","ESCCOR","EVDW2_14",""/
 229       data wname /
 230      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
 231      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
 232      &   "WHPB","WVDWPP","WBOND","WSCCOR","WSCP14",""/
 233       data nprint_ene /18/
 234       data print_order /1,2,3,17,11,12,13,14,4,5,6,7,8,9,10,16,15,18,19,
 235      &   20/
 236 c Dielectric constant of water
 237       data eps_out /80.0d0/
 238       end
 239 c---------------------------------------------------------------------------
 240       subroutine init_int_table
 241       implicit real*8 (a-h,o-z)
 242       include 'DIMENSIONS'
 243       include 'sizesclu.dat'
 244 #ifdef MPI
 245       include 'mpif.h'
 246 #endif
 247 #ifdef MPL
 248       include 'COMMON.INFO'
 249 #endif
 250       include 'COMMON.CHAIN'
 251       include 'COMMON.INTERACT'
 252       include 'COMMON.LOCAL'
 253       include 'COMMON.SBRIDGE'
 254       include 'COMMON.IOUNITS'
 255       logical scheck,lprint
 256 #ifdef MPL
 257       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
 258      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
 259 C... Determine the numbers of start and end SC-SC interaction
 260 C... to deal with by current processor.
 261       lprint=.false.
 262       if (lprint)
 263      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
 264       n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
 265       MyRank=MyID-(MyGroup-1)*fgProcs
 266       call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
 267       if (lprint)
 268      &  write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 269      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
 270      &  ' my_sc_inde',my_sc_inde
 271       ind_sctint=0
 272       iatsc_s=0
 273       iatsc_e=0
 274 #endif
 275       lprint=.false.
 276       do i=1,maxres
 277         nint_gr(i)=0
 278         nscp_gr(i)=0
 279         do j=1,maxint_gr
 280           istart(i,1)=0
 281           iend(i,1)=0
 282           ielstart(i)=0
 283           ielend(i)=0
 284           iscpstart(i,1)=0
 285           iscpend(i,1)=0
 286         enddo
 287       enddo
 288       ind_scint=0
 289       ind_scint_old=0
 290 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
 291 cd   &   (ihpb(i),jhpb(i),i=1,nss)
 292       do i=nnt,nct-1
 293         scheck=.false.
 294         do ii=1,nss
 295           if (ihpb(ii).eq.i+nres) then
 296             scheck=.true.
 297             jj=jhpb(ii)-nres
 298             goto 10
 299           endif
 300         enddo
 301    10   continue
 302 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
 303         if (scheck) then
 304           if (jj.eq.i+1) then
 305 #ifdef MPL
 306             write (iout,*) 'jj=i+1'
 307             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 308      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 309 #else
 310             nint_gr(i)=1
 311             istart(i,1)=i+2
 312             iend(i,1)=nct
 313 #endif
 314           else if (jj.eq.nct) then
 315 #ifdef MPL
 316             write (iout,*) 'jj=nct'
 317             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 318      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 319 #else
 320             nint_gr(i)=1
 321             istart(i,1)=i+1
 322             iend(i,1)=nct-1
 323 #endif
 324           else
 325 #ifdef MPL
 326             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 327      & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 328             ii=nint_gr(i)+1
 329             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 330      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
 331 #else
 332             nint_gr(i)=2
 333             istart(i,1)=i+1
 334             iend(i,1)=jj-1
 335             istart(i,2)=jj+1
 336             iend(i,2)=nct
 337 #endif
 338           endif
 339         else
 340 #ifdef MPL
 341           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 342      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 343 #else
 344           nint_gr(i)=1
 345           istart(i,1)=i+1
 346           iend(i,1)=nct
 347           ind_scint=int_scint+nct-i
 348 #endif
 349         endif
 350 #ifdef MPL
 351         ind_scint_old=ind_scint
 352 #endif
 353       enddo
 354    12 continue
 355 #ifndef MPL
 356       iatsc_s=nnt
 357       iatsc_e=nct-1
 358 #endif
 359 #ifdef MPL
 360       if (lprint) then
 361         write (iout,*) 'Processor',MyID,' Group',MyGroup
 362         write (iout,*) 'iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
 363       endif
 364 #endif
 365       if (lprint) then
 366       write (iout,'(a)') 'Interaction array:'
 367       do i=iatsc_s,iatsc_e
 368         write (iout,'(i3,2(2x,2i3))')
 369      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
 370       enddo
 371       endif
 372       ispp=2
 373 #ifdef MPL
 374 C Now partition the electrostatic-interaction array
 375       npept=nct-nnt
 376       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
 377       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
 378       if (lprint)
 379      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 380      &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
 381      &               ' my_ele_inde',my_ele_inde
 382       iatel_s=0
 383       iatel_e=0
 384       ind_eleint=0
 385       ind_eleint_old=0
 386       do i=nnt,nct-3
 387         ijunk=0
 388         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
 389      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
 390       enddo ! i
 391    13 continue
 392 #else
 393       iatel_s=nnt
 394       iatel_e=nct-3
 395       do i=iatel_s,iatel_e
 396         ielstart(i)=i+2
 397         ielend(i)=nct-1
 398       enddo
 399 #endif
 400       if (lprint) then
 401         write (iout,'(a)') 'Electrostatic interaction array:'
 402         do i=iatel_s,iatel_e
 403           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
 404         enddo
 405       endif ! lprint
 406 c     iscp=3
 407       iscp=2
 408 C Partition the SC-p interaction array
 409 #ifdef MPL
 410       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
 411       call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
 412       if (lprint)
 413      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 414      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
 415      &               ' my_scp_inde',my_scp_inde
 416       iatscp_s=0
 417       iatscp_e=0
 418       ind_scpint=0
 419       ind_scpint_old=0
 420       do i=nnt,nct-1
 421         if (i.lt.nnt+iscp) then
 422 cd        write (iout,*) 'i.le.nnt+iscp'
 423           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 424      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
 425      &      iscpend(i,1),*14)
 426         else if (i.gt.nct-iscp) then
 427 cd        write (iout,*) 'i.gt.nct-iscp'
 428           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 429      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 430      &      iscpend(i,1),*14)
 431         else
 432           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 433      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 434      &      iscpend(i,1),*14)
 435           ii=nscp_gr(i)+1
 436           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 437      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
 438      &      iscpend(i,ii),*14)
 439         endif
 440       enddo ! i
 441    14 continue
 442 #else
 443       iatscp_s=nnt
 444       iatscp_e=nct-1
 445       do i=nnt,nct-1
 446         if (i.lt.nnt+iscp) then
 447           nscp_gr(i)=1
 448           iscpstart(i,1)=i+iscp
 449           iscpend(i,1)=nct
 450         elseif (i.gt.nct-iscp) then
 451           nscp_gr(i)=1
 452           iscpstart(i,1)=nnt
 453           iscpend(i,1)=i-iscp
 454         else
 455           nscp_gr(i)=2
 456           iscpstart(i,1)=nnt
 457           iscpend(i,1)=i-iscp
 458           iscpstart(i,2)=i+iscp
 459           iscpend(i,2)=nct
 460         endif
 461       enddo ! i
 462 #endif
 463       if (lprint) then
 464         write (iout,'(a)') 'SC-p interaction array:'
 465         do i=iatscp_s,iatscp_e
 466           write (iout,'(i3,2(2x,2i3))')
 467      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
 468         enddo
 469       endif ! lprint
 470 C Partition local interactions
 471 #ifdef MPL
 472       call int_bounds(nres-2,loc_start,loc_end)
 473       loc_start=loc_start+1
 474       loc_end=loc_end+1
 475       call int_bounds(nres-2,ithet_start,ithet_end)
 476       ithet_start=ithet_start+2
 477       ithet_end=ithet_end+2
 478       call int_bounds(nct-nnt-2,iphi_start,iphi_end)
 479       iphi_start=iphi_start+nnt+2
 480       iphi_end=iphi_end+nnt+2
 481       call int_bounds(nres-3,itau_start,itau_end)
 482       itau_start=itau_start+3
 483       itau_end=itau_end+3
 484       if (lprint) then
 485         write (iout,*) 'Processor:',MyID,
 486      & ' loc_start',loc_start,' loc_end',loc_end,
 487      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 488      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 489         write (*,*) 'Processor:',MyID,
 490      & ' loc_start',loc_start,' loc_end',loc_end,
 491      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 492      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 493       endif
 494       if (fgprocs.gt.1 .and. MyID.eq.BossID) then
 495         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
 496      & nele_int_tot,' electrostatic and ',nscp_int_tot,
 497      & ' SC-p interactions','were distributed among',fgprocs,
 498      & ' fine-grain processors.'
 499       endif
 500 #else
 501       loc_start=2
 502       loc_end=nres-1
 503       ithet_start=3
 504       ithet_end=nres
 505       iphi_start=nnt+3
 506       iphi_end=nct
 507       itau_start=4
 508       itau_end=nres
 509
 510 #endif
 511       return
 512       end
 513 c---------------------------------------------------------------------------
 514       subroutine int_partition(int_index,lower_index,upper_index,atom,
 515      & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
 516       implicit real*8 (a-h,o-z)
 517       include 'DIMENSIONS'
 518       include 'COMMON.IOUNITS'
 519       integer int_index,lower_index,upper_index,atom,at_start,at_end,
 520      & first_atom,last_atom,int_gr,jat_start,jat_end
 521       logical lprn
 522       lprn=.false.
 523       if (lprn) write (iout,*) 'int_index=',int_index
 524       int_index_old=int_index
 525       int_index=int_index+last_atom-first_atom+1
 526       if (lprn)
 527      &   write (iout,*) 'int_index=',int_index,
 528      &               ' int_index_old',int_index_old,
 529      &               ' lower_index=',lower_index,
 530      &               ' upper_index=',upper_index,
 531      &               ' atom=',atom,' first_atom=',first_atom,
 532      &               ' last_atom=',last_atom
 533       if (int_index.ge.lower_index) then
 534         int_gr=int_gr+1
 535         if (at_start.eq.0) then
 536           at_start=atom
 537           jat_start=first_atom-1+lower_index-int_index_old
 538         else
 539           jat_start=first_atom
 540         endif
 541         if (lprn) write (iout,*) 'jat_start',jat_start
 542         if (int_index.ge.upper_index) then
 543           at_end=atom
 544           jat_end=first_atom-1+upper_index-int_index_old
 545           return1
 546         else
 547           jat_end=last_atom
 548         endif
 549         if (lprn) write (iout,*) 'jat_end',jat_end
 550       endif
 551       return
 552       end
 553 c------------------------------------------------------------------------------
 554       subroutine hpb_partition
 555       implicit real*8 (a-h,o-z)
 556       include 'DIMENSIONS'
 557       include 'COMMON.SBRIDGE'
 558       include 'COMMON.IOUNITS'
 559 #ifdef MPL
 560       include 'COMMON.INFO'
 561       call int_bounds(nhpb,link_start,link_end)
 562 #else
 563       link_start=1
 564       link_end=nhpb
 565 #endif
 566 cd    write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 567 cd   &  ' nhpb',nhpb,' link_start=',link_start,
 568 cd   &  ' link_end',link_end
 569       return
 570       end