source/cluster/wham/src/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          do ichir1=-1,1
 105           do ichir2=-1,1
 106           athet(j,i,ichir1,ichir2)=0.0D0
 107           bthet(j,i,ichir1,ichir2)=0.0D0
 108           enddo
 109          enddo
 110         enddo
 111         do j=0,3
 112           polthet(j,i)=0.0D0
 113         enddo
 114         do j=1,3
 115           gthet(j,i)=0.0D0
 116         enddo
 117         theta0(i)=0.0D0
 118         sig0(i)=0.0D0
 119         sigc0(i)=0.0D0
 120         do j=1,maxlob
 121           bsc(j,i)=0.0D0
 122           do k=1,3
 123             censc(k,j,i)=0.0D0
 124           enddo
 125           do k=1,3
 126             do l=1,3
 127               gaussc(l,k,j,i)=0.0D0
 128             enddo
 129           enddo
 130           nlob(i)=0
 131         enddo
 132       enddo
 133       nlob(ntyp1)=0
 134       dsc(ntyp1)=0.0D0
 135       do i=-maxtor,maxtor
 136         itortyp(i)=0
 137        do iblock=1,2
 138         do j=-maxtor,maxtor
 139           do k=1,maxterm
 140             v1(k,j,i,iblock)=0.0D0
 141             v2(k,j,i,iblock)=0.0D0
 142            enddo
 143          enddo
 144         enddo
 145       enddo
 146       do iblock=1,2
 147        do i=-maxtor,maxtor
 148         do j=-maxtor,maxtor
 149          do k=-maxtor,maxtor
 150           do l=1,maxtermd_1
 151             v1c(1,l,i,j,k,iblock)=0.0D0
 152             v1s(1,l,i,j,k,iblock)=0.0D0
 153             v1c(2,l,i,j,k,iblock)=0.0D0
 154             v1s(2,l,i,j,k,iblock)=0.0D0
 155           enddo !l
 156           do l=1,maxtermd_2
 157            do m=1,maxtermd_2
 158             v2c(m,l,i,j,k,iblock)=0.0D0
 159             v2s(m,l,i,j,k,iblock)=0.0D0
 160            enddo !m
 161           enddo !l
 162         enddo !k
 163        enddo !j
 164       enddo !i
 165       enddo !iblock
 166       do i=1,maxres
 167         itype(i)=0
 168         itel(i)=0
 169       enddo
 170 C Initialize the bridge arrays
 171       ns=0
 172       nss=0
 173       nhpb=0
 174       do i=1,maxss
 175         iss(i)=0
 176       enddo
 177       do i=1,maxss
 178         dhpb(i)=0.0D0
 179       enddo
 180       do i=1,maxss
 181         ihpb(i)=0
 182         jhpb(i)=0
 183       enddo
 184       do i=1,maxres
 185         dyn_ss_mask(i)=.false.
 186       enddo
 187 C
 188 C Initialize timing.
 189 C
 190       call set_timers
 191 C
 192 C Initialize variables used in minimization.
 193 C
 194 c     maxfun=5000
 195 c     maxit=2000
 196       maxfun=500
 197       maxit=200
 198       tolf=1.0D-2
 199       rtolf=5.0D-4
 200 C
 201 C Initialize the variables responsible for the mode of gradient storage.
 202 C
 203       nfl=0
 204       icg=1
 205       do i=1,14
 206         do j=1,14
 207           if (print_order(i).eq.j) then
 208             iw(print_order(i))=j
 209             goto 1121
 210           endif
 211         enddo
 212 1121    continue
 213       enddo
 214       calc_grad=.false.
 215 C Set timers and counters for the respective routines
 216       t_func = 0.0d0
 217       t_grad = 0.0d0
 218       t_fhel = 0.0d0
 219       t_fbet = 0.0d0
 220       t_ghel = 0.0d0
 221       t_gbet = 0.0d0
 222       t_viol = 0.0d0
 223       t_gviol = 0.0d0
 224       n_func = 0
 225       n_grad = 0
 226       n_fhel = 0
 227       n_fbet = 0
 228       n_ghel = 0
 229       n_gbet = 0
 230       n_viol = 0
 231       n_gviol = 0
 232       n_map = 0
 233 #ifndef SPLITELE
 234       nprint_ene=nprint_ene-1
 235 #endif
 236       return
 237       end
 238 c-------------------------------------------------------------------------
 239       block data nazwy
 240       implicit real*8 (a-h,o-z)
 241       include 'DIMENSIONS'
 242       include 'sizesclu.dat'
 243       include 'COMMON.NAMES'
 244       include 'COMMON.FFIELD'
 245       data restyp /
 246      &'DD' ,'DPR','DLY','DAR','DHI','DAS','DGL','DSG','DGN','DSN','DTH',
 247      &'DYY','DAL','DTY','DTR','DVA','DLE','DIL','DPN','MED','DCY','ZER',
 248      &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
 249      &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','D'/
 250       data onelet /
 251      &'z','p','k','r','h','d','e','n','q','s','t','g',
 252      &'a','y','w','v','l','i','f','m','c','x',
 253      &'C','M','F','I','L','V','W','Y','A','G','T',
 254      &'S','Q','N','E','D','H','R','K','P','X'/
 255       data potname /'LJ','LJK','BP','GB','GBV'/
 256       data ename /
 257      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
 258      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
 259      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP",
 260      &   "ESTR","ESCCOR","EVDW2_14",""/
 261       data wname /
 262      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
 263      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
 264      &   "WHPB","WVDWPP","WBOND","WSCCOR","WSCP14",""/
 265       data nprint_ene /18/
 266       data print_order /1,2,3,17,11,12,13,14,4,5,6,7,8,9,10,16,15,18,19,
 267      &   20/
 268       end
 269 c---------------------------------------------------------------------------
 270       subroutine init_int_table
 271       implicit real*8 (a-h,o-z)
 272       include 'DIMENSIONS'
 273       include 'sizesclu.dat'
 274 #ifdef MPI
 275       include 'mpif.h'
 276 #endif
 277 #ifdef MPL
 278       include 'COMMON.INFO'
 279 #endif
 280       include 'COMMON.CHAIN'
 281       include 'COMMON.INTERACT'
 282       include 'COMMON.LOCAL'
 283       include 'COMMON.SBRIDGE'
 284       include 'COMMON.IOUNITS'
 285       logical scheck,lprint
 286 #ifdef MPL
 287       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
 288      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
 289 C... Determine the numbers of start and end SC-SC interaction
 290 C... to deal with by current processor.
 291       lprint=.false.
 292       if (lprint)
 293      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
 294       n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
 295       MyRank=MyID-(MyGroup-1)*fgProcs
 296       call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
 297       if (lprint)
 298      &  write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 299      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
 300      &  ' my_sc_inde',my_sc_inde
 301       ind_sctint=0
 302       iatsc_s=0
 303       iatsc_e=0
 304 #endif
 305       lprint=.false.
 306       do i=1,maxres
 307         nint_gr(i)=0
 308         nscp_gr(i)=0
 309         do j=1,maxint_gr
 310           istart(i,1)=0
 311           iend(i,1)=0
 312           ielstart(i)=0
 313           ielend(i)=0
 314           iscpstart(i,1)=0
 315           iscpend(i,1)=0
 316         enddo
 317       enddo
 318       ind_scint=0
 319       ind_scint_old=0
 320 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
 321 cd   &   (ihpb(i),jhpb(i),i=1,nss)
 322       do i=nnt,nct-1
 323         scheck=.false.
 324         do ii=1,nss
 325           if (ihpb(ii).eq.i+nres) then
 326             scheck=.true.
 327             if (dyn_ss) go to 10
 328             jj=jhpb(ii)-nres
 329             goto 10
 330           endif
 331         enddo
 332    10   continue
 333 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
 334         if (scheck) then
 335           if (jj.eq.i+1) then
 336 #ifdef MPL
 337             write (iout,*) 'jj=i+1'
 338             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 339      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 340 #else
 341             nint_gr(i)=1
 342             istart(i,1)=i+2
 343             iend(i,1)=nct
 344 #endif
 345           else if (jj.eq.nct) then
 346 #ifdef MPL
 347             write (iout,*) 'jj=nct'
 348             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 349      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 350 #else
 351             nint_gr(i)=1
 352             istart(i,1)=i+1
 353             iend(i,1)=nct-1
 354 #endif
 355           else
 356 #ifdef MPL
 357             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 358      & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 359             ii=nint_gr(i)+1
 360             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 361      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
 362 #else
 363             nint_gr(i)=2
 364             istart(i,1)=i+1
 365             iend(i,1)=jj-1
 366             istart(i,2)=jj+1
 367             iend(i,2)=nct
 368 #endif
 369           endif
 370         else
 371 #ifdef MPL
 372           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 373      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 374 #else
 375           nint_gr(i)=1
 376           istart(i,1)=i+1
 377           iend(i,1)=nct
 378           ind_scint=int_scint+nct-i
 379 #endif
 380         endif
 381 #ifdef MPL
 382         ind_scint_old=ind_scint
 383 #endif
 384       enddo
 385    12 continue
 386 #ifndef MPL
 387       iatsc_s=nnt
 388       iatsc_e=nct-1
 389 #endif
 390 #ifdef MPL
 391       if (lprint) then
 392         write (iout,*) 'Processor',MyID,' Group',MyGroup
 393         write (iout,*) 'iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
 394       endif
 395 #endif
 396       if (lprint) then
 397       write (iout,'(a)') 'Interaction array:'
 398       do i=iatsc_s,iatsc_e
 399         write (iout,'(i3,2(2x,2i3))')
 400      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
 401       enddo
 402       endif
 403       ispp=2
 404 #ifdef MPL
 405 C Now partition the electrostatic-interaction array
 406       npept=nct-nnt
 407       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
 408       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
 409       if (lprint)
 410      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 411      &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
 412      &               ' my_ele_inde',my_ele_inde
 413       iatel_s=0
 414       iatel_e=0
 415       ind_eleint=0
 416       ind_eleint_old=0
 417       do i=nnt,nct-3
 418         ijunk=0
 419         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
 420      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
 421       enddo ! i
 422    13 continue
 423 #else
 424       iatel_s=nnt
 425       iatel_e=nct-3
 426       do i=iatel_s,iatel_e
 427         ielstart(i)=i+2
 428         ielend(i)=nct-1
 429       enddo
 430 #endif
 431       if (lprint) then
 432         write (iout,'(a)') 'Electrostatic interaction array:'
 433         do i=iatel_s,iatel_e
 434           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
 435         enddo
 436       endif ! lprint
 437 c     iscp=3
 438       iscp=2
 439 C Partition the SC-p interaction array
 440 #ifdef MPL
 441       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
 442       call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
 443       if (lprint)
 444      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 445      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
 446      &               ' my_scp_inde',my_scp_inde
 447       iatscp_s=0
 448       iatscp_e=0
 449       ind_scpint=0
 450       ind_scpint_old=0
 451       do i=nnt,nct-1
 452         if (i.lt.nnt+iscp) then
 453 cd        write (iout,*) 'i.le.nnt+iscp'
 454           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 455      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
 456      &      iscpend(i,1),*14)
 457         else if (i.gt.nct-iscp) then
 458 cd        write (iout,*) 'i.gt.nct-iscp'
 459           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 460      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 461      &      iscpend(i,1),*14)
 462         else
 463           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 464      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 465      &      iscpend(i,1),*14)
 466           ii=nscp_gr(i)+1
 467           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 468      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
 469      &      iscpend(i,ii),*14)
 470         endif
 471       enddo ! i
 472    14 continue
 473 #else
 474       iatscp_s=nnt
 475       iatscp_e=nct-1
 476       do i=nnt,nct-1
 477         if (i.lt.nnt+iscp) then
 478           nscp_gr(i)=1
 479           iscpstart(i,1)=i+iscp
 480           iscpend(i,1)=nct
 481         elseif (i.gt.nct-iscp) then
 482           nscp_gr(i)=1
 483           iscpstart(i,1)=nnt
 484           iscpend(i,1)=i-iscp
 485         else
 486           nscp_gr(i)=2
 487           iscpstart(i,1)=nnt
 488           iscpend(i,1)=i-iscp
 489           iscpstart(i,2)=i+iscp
 490           iscpend(i,2)=nct
 491         endif
 492       enddo ! i
 493 #endif
 494       if (lprint) then
 495         write (iout,'(a)') 'SC-p interaction array:'
 496         do i=iatscp_s,iatscp_e
 497           write (iout,'(i3,2(2x,2i3))')
 498      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
 499         enddo
 500       endif ! lprint
 501 C Partition local interactions
 502 #ifdef MPL
 503       call int_bounds(nres-2,loc_start,loc_end)
 504       loc_start=loc_start+1
 505       loc_end=loc_end+1
 506       call int_bounds(nres-2,ithet_start,ithet_end)
 507       ithet_start=ithet_start+2
 508       ithet_end=ithet_end+2
 509       call int_bounds(nct-nnt-2,iphi_start,iphi_end)
 510       iphi_start=iphi_start+nnt+2
 511       iphi_end=iphi_end+nnt+2
 512       call int_bounds(nres-3,itau_start,itau_end)
 513       itau_start=itau_start+3
 514       itau_end=itau_end+3
 515       if (lprint) then
 516         write (iout,*) 'Processor:',MyID,
 517      & ' loc_start',loc_start,' loc_end',loc_end,
 518      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 519      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 520         write (*,*) 'Processor:',MyID,
 521      & ' loc_start',loc_start,' loc_end',loc_end,
 522      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 523      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 524       endif
 525       if (fgprocs.gt.1 .and. MyID.eq.BossID) then
 526         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
 527      & nele_int_tot,' electrostatic and ',nscp_int_tot,
 528      & ' SC-p interactions','were distributed among',fgprocs,
 529      & ' fine-grain processors.'
 530       endif
 531 #else
 532       loc_start=2
 533       loc_end=nres-1
 534       ithet_start=3
 535       ithet_end=nres
 536       iphi_start=nnt+3
 537       iphi_end=nct
 538       itau_start=4
 539       itau_end=nres
 540
 541 #endif
 542       return
 543       end
 544 c---------------------------------------------------------------------------
 545       subroutine int_partition(int_index,lower_index,upper_index,atom,
 546      & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
 547       implicit real*8 (a-h,o-z)
 548       include 'DIMENSIONS'
 549       include 'COMMON.IOUNITS'
 550       integer int_index,lower_index,upper_index,atom,at_start,at_end,
 551      & first_atom,last_atom,int_gr,jat_start,jat_end
 552       logical lprn
 553       lprn=.false.
 554       if (lprn) write (iout,*) 'int_index=',int_index
 555       int_index_old=int_index
 556       int_index=int_index+last_atom-first_atom+1
 557       if (lprn)
 558      &   write (iout,*) 'int_index=',int_index,
 559      &               ' int_index_old',int_index_old,
 560      &               ' lower_index=',lower_index,
 561      &               ' upper_index=',upper_index,
 562      &               ' atom=',atom,' first_atom=',first_atom,
 563      &               ' last_atom=',last_atom
 564       if (int_index.ge.lower_index) then
 565         int_gr=int_gr+1
 566         if (at_start.eq.0) then
 567           at_start=atom
 568           jat_start=first_atom-1+lower_index-int_index_old
 569         else
 570           jat_start=first_atom
 571         endif
 572         if (lprn) write (iout,*) 'jat_start',jat_start
 573         if (int_index.ge.upper_index) then
 574           at_end=atom
 575           jat_end=first_atom-1+upper_index-int_index_old
 576           return1
 577         else
 578           jat_end=last_atom
 579         endif
 580         if (lprn) write (iout,*) 'jat_end',jat_end
 581       endif
 582       return
 583       end
 584 c------------------------------------------------------------------------------
 585       subroutine hpb_partition
 586       implicit real*8 (a-h,o-z)
 587       include 'DIMENSIONS'
 588       include 'COMMON.SBRIDGE'
 589       include 'COMMON.IOUNITS'
 590 #ifdef MPL
 591       include 'COMMON.INFO'
 592       call int_bounds(nhpb,link_start,link_end)
 593 #else
 594       link_start=1
 595       link_end=nhpb
 596 #endif
 597 cd    write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 598 cd   &  ' nhpb',nhpb,' link_start=',link_start,
 599 cd   &  ' link_end',link_end
 600       return
 601       end