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          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 C
 185 C Initialize timing.
 186 C
 187       call set_timers
 188 C
 189 C Initialize variables used in minimization.
 190 C
 191 c     maxfun=5000
 192 c     maxit=2000
 193       maxfun=500
 194       maxit=200
 195       tolf=1.0D-2
 196       rtolf=5.0D-4
 197 C
 198 C Initialize the variables responsible for the mode of gradient storage.
 199 C
 200       nfl=0
 201       icg=1
 202       do i=1,14
 203         do j=1,14
 204           if (print_order(i).eq.j) then
 205             iw(print_order(i))=j
 206             goto 1121
 207           endif
 208         enddo
 209 1121    continue
 210       enddo
 211       calc_grad=.false.
 212 C Set timers and counters for the respective routines
 213       t_func = 0.0d0
 214       t_grad = 0.0d0
 215       t_fhel = 0.0d0
 216       t_fbet = 0.0d0
 217       t_ghel = 0.0d0
 218       t_gbet = 0.0d0
 219       t_viol = 0.0d0
 220       t_gviol = 0.0d0
 221       n_func = 0
 222       n_grad = 0
 223       n_fhel = 0
 224       n_fbet = 0
 225       n_ghel = 0
 226       n_gbet = 0
 227       n_viol = 0
 228       n_gviol = 0
 229       n_map = 0
 230 #ifndef SPLITELE
 231       nprint_ene=nprint_ene-1
 232 #endif
 233       return
 234       end
 235 c-------------------------------------------------------------------------
 236       block data nazwy
 237       implicit real*8 (a-h,o-z)
 238       include 'DIMENSIONS'
 239       include 'sizesclu.dat'
 240       include 'COMMON.NAMES'
 241       include 'COMMON.FFIELD'
 242       data restyp /
 243      &'DD','DAU','DAI','DDB','DSM','DPR','DLY','DAR','DHI','DAS','DGL',
 244      & 'DSG','DGN','DSN','DTH',
 245      &'DYY','DAL','DTY','DTR','DVA','DLE','DIL','DPN','MED','DCY','ZER',
 246      &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
 247      &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','SME','DBZ',
 248      &'AIB','ABU','D'/
 249       data onelet /
 250      &'z','z','z','z','z','p','k','r','h','d','e','n','q','s','t','g',
 251      &'a','y','w','v','l','i','f','m','c','x',
 252      &'C','M','F','I','L','V','W','Y','A','G','T',
 253      &'S','Q','N','E','D','H','R','K','P','z','z','z','z','X'/
 254       data potname /'LJ','LJK','BP','GB','GBV'/
 255       data ename /
 256      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
 257      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
 258      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP",
 259      &   "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T","ELIPTRAN",
 260      &   "EAFM","ETHETC"/
 261       data wname /
 262      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
 263      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
 264      &   "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC",
 265      &   "WLIPTRAN","WAFM","WTHETC"/
 266       data nprint_ene /22/
 267       data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,
 268      &  16,15,17,20,21,24,22,23/
 269       end
 270 c---------------------------------------------------------------------------
 271       subroutine init_int_table
 272       implicit real*8 (a-h,o-z)
 273       include 'DIMENSIONS'
 274       include 'sizesclu.dat'
 275 #ifdef MPI
 276       include 'mpif.h'
 277 #endif
 278 #ifdef MPL
 279       include 'COMMON.INFO'
 280 #endif
 281       include 'COMMON.CHAIN'
 282       include 'COMMON.INTERACT'
 283       include 'COMMON.LOCAL'
 284       include 'COMMON.SBRIDGE'
 285       include 'COMMON.IOUNITS'
 286       logical scheck,lprint
 287 #ifdef MPL
 288       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
 289      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
 290 C... Determine the numbers of start and end SC-SC interaction
 291 C... to deal with by current processor.
 292       lprint=.false.
 293       if (lprint)
 294      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
 295       n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
 296       MyRank=MyID-(MyGroup-1)*fgProcs
 297       call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
 298       if (lprint)
 299      &  write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 300      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
 301      &  ' my_sc_inde',my_sc_inde
 302       ind_sctint=0
 303       iatsc_s=0
 304       iatsc_e=0
 305 #endif
 306       lprint=.false.
 307       do i=1,maxres
 308         nint_gr(i)=0
 309         nscp_gr(i)=0
 310         do j=1,maxint_gr
 311           istart(i,1)=0
 312           iend(i,1)=0
 313           ielstart(i)=0
 314           ielend(i)=0
 315           iscpstart(i,1)=0
 316           iscpend(i,1)=0
 317         enddo
 318       enddo
 319       ind_scint=0
 320       ind_scint_old=0
 321 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
 322 cd   &   (ihpb(i),jhpb(i),i=1,nss)
 323       do i=nnt,nct-1
 324         scheck=.false.
 325         if (dyn_ss) goto 10
 326         do ii=1,nss
 327           if (ihpb(ii).eq.i+nres) then
 328             scheck=.true.
 329             jj=jhpb(ii)-nres
 330             goto 10
 331           endif
 332         enddo
 333    10   continue
 334 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
 335         if (scheck) then
 336           if (jj.eq.i+1) then
 337 #ifdef MPL
 338             write (iout,*) 'jj=i+1'
 339             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 340      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 341 #else
 342             nint_gr(i)=1
 343             istart(i,1)=i+2
 344             iend(i,1)=nct
 345 #endif
 346           else if (jj.eq.nct) then
 347 #ifdef MPL
 348             write (iout,*) 'jj=nct'
 349             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 350      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 351 #else
 352             nint_gr(i)=1
 353             istart(i,1)=i+1
 354             iend(i,1)=nct-1
 355 #endif
 356           else
 357 #ifdef MPL
 358             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 359      & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 360             ii=nint_gr(i)+1
 361             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 362      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
 363 #else
 364             nint_gr(i)=2
 365             istart(i,1)=i+1
 366             iend(i,1)=jj-1
 367             istart(i,2)=jj+1
 368             iend(i,2)=nct
 369 #endif
 370           endif
 371         else
 372 #ifdef MPL
 373           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 374      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 375 #else
 376           nint_gr(i)=1
 377           istart(i,1)=i+1
 378           iend(i,1)=nct
 379           ind_scint=ind_scint+nct-i
 380 #endif
 381         endif
 382 #ifdef MPL
 383         ind_scint_old=ind_scint
 384 #endif
 385       enddo
 386    12 continue
 387 #ifndef MPL
 388       iatsc_s=nnt
 389       iatsc_e=nct-1
 390 #endif
 391 #ifdef MPL
 392       if (lprint) then
 393         write (iout,*) 'Processor',MyID,' Group',MyGroup
 394         write (iout,*) 'iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
 395       endif
 396 #endif
 397       if (lprint) then
 398       write (iout,'(a)') 'Interaction array:'
 399       do i=iatsc_s,iatsc_e
 400         write (iout,'(i3,2(2x,2i3))')
 401      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
 402       enddo
 403       endif
 404       ispp=2
 405 #ifdef MPL
 406 C Now partition the electrostatic-interaction array
 407       npept=nct-nnt
 408       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
 409       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
 410       if (lprint)
 411      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 412      &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
 413      &               ' my_ele_inde',my_ele_inde
 414       iatel_s=0
 415       iatel_e=0
 416       ind_eleint=0
 417       ind_eleint_old=0
 418       do i=nnt,nct-3
 419         ijunk=0
 420         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
 421      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
 422       enddo ! i
 423    13 continue
 424 #else
 425       iatel_s=nnt
 426       iatel_e=nct-3
 427       do i=iatel_s,iatel_e
 428         ielstart(i)=i+2
 429         ielend(i)=nct-1
 430       enddo
 431 #endif
 432       if (lprint) then
 433         write (iout,'(a)') 'Electrostatic interaction array:'
 434         do i=iatel_s,iatel_e
 435           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
 436         enddo
 437       endif ! lprint
 438 c     iscp=3
 439       iscp=2
 440 C Partition the SC-p interaction array
 441 #ifdef MPL
 442       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
 443       call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
 444       if (lprint)
 445      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 446      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
 447      &               ' my_scp_inde',my_scp_inde
 448       iatscp_s=0
 449       iatscp_e=0
 450       ind_scpint=0
 451       ind_scpint_old=0
 452       do i=nnt,nct-1
 453         if (i.lt.nnt+iscp) then
 454 cd        write (iout,*) 'i.le.nnt+iscp'
 455           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 456      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
 457      &      iscpend(i,1),*14)
 458         else if (i.gt.nct-iscp) then
 459 cd        write (iout,*) 'i.gt.nct-iscp'
 460           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 461      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 462      &      iscpend(i,1),*14)
 463         else
 464           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 465      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 466      &      iscpend(i,1),*14)
 467           ii=nscp_gr(i)+1
 468           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 469      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
 470      &      iscpend(i,ii),*14)
 471         endif
 472       enddo ! i
 473    14 continue
 474 #else
 475       iatscp_s=nnt
 476       iatscp_e=nct-1
 477       do i=nnt,nct-1
 478         if (i.lt.nnt+iscp) then
 479           nscp_gr(i)=1
 480           iscpstart(i,1)=i+iscp
 481           iscpend(i,1)=nct
 482         elseif (i.gt.nct-iscp) then
 483           nscp_gr(i)=1
 484           iscpstart(i,1)=nnt
 485           iscpend(i,1)=i-iscp
 486         else
 487           nscp_gr(i)=2
 488           iscpstart(i,1)=nnt
 489           iscpend(i,1)=i-iscp
 490           iscpstart(i,2)=i+iscp
 491           iscpend(i,2)=nct
 492         endif
 493       enddo ! i
 494 #endif
 495       if (lprint) then
 496         write (iout,'(a)') 'SC-p interaction array:'
 497         do i=iatscp_s,iatscp_e
 498           write (iout,'(i3,2(2x,2i3))')
 499      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
 500         enddo
 501       endif ! lprint
 502 C Partition local interactions
 503 #ifdef MPL
 504       call int_bounds(nres-2,loc_start,loc_end)
 505       loc_start=loc_start+1
 506       loc_end=loc_end+1
 507       call int_bounds(nres-2,ithet_start,ithet_end)
 508       ithet_start=ithet_start+2
 509       ithet_end=ithet_end+2
 510       call int_bounds(nct-nnt-2,iphi_start,iphi_end)
 511       iphi_start=iphi_start+nnt+2
 512       iphi_end=iphi_end+nnt+2
 513       call int_bounds(nres-3,itau_start,itau_end)
 514       itau_start=itau_start+3
 515       itau_end=itau_end+3
 516       if (lprint) then
 517         write (iout,*) 'Processor:',MyID,
 518      & ' loc_start',loc_start,' loc_end',loc_end,
 519      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 520      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 521         write (*,*) 'Processor:',MyID,
 522      & ' loc_start',loc_start,' loc_end',loc_end,
 523      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 524      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 525       endif
 526       if (fgprocs.gt.1 .and. MyID.eq.BossID) then
 527         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
 528      & nele_int_tot,' electrostatic and ',nscp_int_tot,
 529      & ' SC-p interactions','were distributed among',fgprocs,
 530      & ' fine-grain processors.'
 531       endif
 532 #else
 533       loc_start=2
 534       loc_end=nres-1
 535       ithet_start=3
 536       ithet_end=nres
 537       iphi_start=nnt+3
 538       iphi_end=nct
 539       itau_start=4
 540       itau_end=nres
 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