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       iliptran=60
  62 C
  63 C Set default weights of the energy terms.
  64 C
  65       wlong=1.0D0
  66       welec=1.0D0
  67       wtor =1.0D0
  68       wang =1.0D0
  69       wscloc=1.0D0
  70       wstrain=1.0D0
  71 C
  72 C Zero out tables.
  73 C
  74       ndih_constr=0
  75       do i=1,maxres2
  76         do j=1,3
  77           c(j,i)=0.0D0
  78           dc(j,i)=0.0D0
  79         enddo
  80       enddo
  81       do i=1,maxres
  82         do j=1,3
  83           xloc(j,i)=0.0D0
  84         enddo
  85       enddo
  86       do i=1,ntyp
  87         do j=1,ntyp
  88           aa_aq(i,j)=0.0D0
  89           bb_aq(i,j)=0.0D0
  90           aa_lip(i,j)=0.0D0
  91           bb_lip(i,j)=0.0D0
  92           augm(i,j)=0.0D0
  93           sigma(i,j)=0.0D0
  94           r0(i,j)=0.0D0
  95           chi(i,j)=0.0D0
  96         enddo
  97         do j=1,2
  98           bad(i,j)=0.0D0
  99         enddo
 100         chip(i)=0.0D0
 101         alp(i)=0.0D0
 102         sigma0(i)=0.0D0
 103         sigii(i)=0.0D0
 104         rr0(i)=0.0D0
 105         a0thet(i)=0.0D0
 106         do j=1,2
 107          do ichir1=-1,1
 108           do ichir2=-1,1
 109           athet(j,i,ichir1,ichir2)=0.0D0
 110           bthet(j,i,ichir1,ichir2)=0.0D0
 111           enddo
 112          enddo
 113         enddo
 114         do j=0,3
 115           polthet(j,i)=0.0D0
 116         enddo
 117         do j=1,3
 118           gthet(j,i)=0.0D0
 119         enddo
 120         theta0(i)=0.0D0
 121         sig0(i)=0.0D0
 122         sigc0(i)=0.0D0
 123         do j=1,maxlob
 124           bsc(j,i)=0.0D0
 125           do k=1,3
 126             censc(k,j,i)=0.0D0
 127           enddo
 128           do k=1,3
 129             do l=1,3
 130               gaussc(l,k,j,i)=0.0D0
 131             enddo
 132           enddo
 133           nlob(i)=0
 134         enddo
 135       enddo
 136       nlob(ntyp1)=0
 137       dsc(ntyp1)=0.0D0
 138       do i=-maxtor,maxtor
 139         itortyp(i)=0
 140        do iblock=1,2
 141         do j=-maxtor,maxtor
 142           do k=1,maxterm
 143             v1(k,j,i,iblock)=0.0D0
 144             v2(k,j,i,iblock)=0.0D0
 145            enddo
 146          enddo
 147        enddo
 148       enddo
 149       do iblock=1,2
 150        do i=-maxtor,maxtor
 151         do j=-maxtor,maxtor
 152          do k=-maxtor,maxtor
 153           do l=1,maxtermd_1
 154             v1c(1,l,i,j,k,iblock)=0.0D0
 155             v1s(1,l,i,j,k,iblock)=0.0D0
 156             v1c(2,l,i,j,k,iblock)=0.0D0
 157             v1s(2,l,i,j,k,iblock)=0.0D0
 158           enddo !l
 159           do l=1,maxtermd_2
 160            do m=1,maxtermd_2
 161             v2c(m,l,i,j,k,iblock)=0.0D0
 162             v2s(m,l,i,j,k,iblock)=0.0D0
 163            enddo !m
 164           enddo !l
 165         enddo !k
 166        enddo !j
 167       enddo !i
 168       enddo !iblock
 169       do i=1,maxres
 170         itype(i)=0
 171         itel(i)=0
 172       enddo
 173 C Initialize the bridge arrays
 174       ns=0
 175       nss=0
 176       nhpb=0
 177       do i=1,maxss
 178         iss(i)=0
 179       enddo
 180       do i=1,maxss
 181         dhpb(i)=0.0D0
 182       enddo
 183       do i=1,maxss
 184         ihpb(i)=0
 185         jhpb(i)=0
 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','DAU','DAI','DDB','DSM','DPR','DLY','DAR','DHI','DAS','DGL',
 247      & 'DSG','DGN','DSN','DTH',
 248      &'DYY','DAL','DTY','DTR','DVA','DLE','DIL','DPN','MED','DCY','ZER',
 249      &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
 250      &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','SME','DBZ',
 251      &'AIB','ABU','D'/
 252       data onelet /
 253      &'z','z','z','z','z','p','k','r','h','d','e','n','q','s','t','g',
 254      &'a','y','w','v','l','i','f','m','c','x',
 255      &'C','M','F','I','L','V','W','Y','A','G','T',
 256      &'S','Q','N','E','D','H','R','K','P','z','z','z','z','X'/
 257       data potname /'LJ','LJK','BP','GB','GBV'/
 258       data ename /
 259      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
 260      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
 261      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP",
 262      &   "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T","ELIPTRAN",
 263      &   "EAFM","ETHETC"/
 264       data wname /
 265      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
 266      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
 267      &   "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC",
 268      &   "WLIPTRAN","WAFM","WTHETC"/
 269       data nprint_ene /22/
 270       data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,
 271      &  16,15,17,20,21,24,22,23/
 272       end
 273 c---------------------------------------------------------------------------
 274       subroutine init_int_table
 275       implicit real*8 (a-h,o-z)
 276       include 'DIMENSIONS'
 277       include 'sizesclu.dat'
 278 #ifdef MPI
 279       include 'mpif.h'
 280 #endif
 281 #ifdef MPL
 282       include 'COMMON.INFO'
 283 #endif
 284       include 'COMMON.CHAIN'
 285       include 'COMMON.INTERACT'
 286       include 'COMMON.LOCAL'
 287       include 'COMMON.SBRIDGE'
 288       include 'COMMON.IOUNITS'
 289       logical scheck,lprint
 290 #ifdef MPL
 291       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
 292      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
 293 C... Determine the numbers of start and end SC-SC interaction
 294 C... to deal with by current processor.
 295       lprint=.false.
 296       if (lprint)
 297      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
 298       n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
 299       MyRank=MyID-(MyGroup-1)*fgProcs
 300       call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
 301       if (lprint)
 302      &  write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 303      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
 304      &  ' my_sc_inde',my_sc_inde
 305       ind_sctint=0
 306       iatsc_s=0
 307       iatsc_e=0
 308 #endif
 309       lprint=.false.
 310       do i=1,maxres
 311         nint_gr(i)=0
 312         nscp_gr(i)=0
 313         do j=1,maxint_gr
 314           istart(i,1)=0
 315           iend(i,1)=0
 316           ielstart(i)=0
 317           ielend(i)=0
 318           iscpstart(i,1)=0
 319           iscpend(i,1)=0
 320         enddo
 321       enddo
 322       ind_scint=0
 323       ind_scint_old=0
 324 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
 325 cd   &   (ihpb(i),jhpb(i),i=1,nss)
 326       do i=nnt,nct-1
 327         scheck=.false.
 328         if (dyn_ss) goto 10
 329         do ii=1,nss
 330           if (ihpb(ii).eq.i+nres) then
 331             scheck=.true.
 332             jj=jhpb(ii)-nres
 333             goto 10
 334           endif
 335         enddo
 336    10   continue
 337 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
 338         if (scheck) then
 339           if (jj.eq.i+1) then
 340 #ifdef MPL
 341             write (iout,*) 'jj=i+1'
 342             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 343      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 344 #else
 345             nint_gr(i)=1
 346             istart(i,1)=i+2
 347             iend(i,1)=nct
 348 #endif
 349           else if (jj.eq.nct) then
 350 #ifdef MPL
 351             write (iout,*) 'jj=nct'
 352             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 353      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 354 #else
 355             nint_gr(i)=1
 356             istart(i,1)=i+1
 357             iend(i,1)=nct-1
 358 #endif
 359           else
 360 #ifdef MPL
 361             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 362      & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 363             ii=nint_gr(i)+1
 364             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 365      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
 366 #else
 367             nint_gr(i)=2
 368             istart(i,1)=i+1
 369             iend(i,1)=jj-1
 370             istart(i,2)=jj+1
 371             iend(i,2)=nct
 372 #endif
 373           endif
 374         else
 375 #ifdef MPL
 376           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 377      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 378 #else
 379           nint_gr(i)=1
 380           istart(i,1)=i+1
 381           iend(i,1)=nct
 382           ind_scint=ind_scint+nct-i
 383 #endif
 384         endif
 385 #ifdef MPL
 386         ind_scint_old=ind_scint
 387 #endif
 388       enddo
 389    12 continue
 390 #ifndef MPL
 391       iatsc_s=nnt
 392       iatsc_e=nct-1
 393 #endif
 394 #ifdef MPL
 395       if (lprint) then
 396         write (iout,*) 'Processor',MyID,' Group',MyGroup
 397         write (iout,*) 'iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
 398       endif
 399 #endif
 400       if (lprint) then
 401       write (iout,'(a)') 'Interaction array:'
 402       do i=iatsc_s,iatsc_e
 403         write (iout,'(i3,2(2x,2i3))')
 404      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
 405       enddo
 406       endif
 407       ispp=2
 408 #ifdef MPL
 409 C Now partition the electrostatic-interaction array
 410       npept=nct-nnt
 411       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
 412       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
 413       if (lprint)
 414      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 415      &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
 416      &               ' my_ele_inde',my_ele_inde
 417       iatel_s=0
 418       iatel_e=0
 419       ind_eleint=0
 420       ind_eleint_old=0
 421       do i=nnt,nct-3
 422         ijunk=0
 423         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
 424      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
 425       enddo ! i
 426    13 continue
 427 #else
 428       iatel_s=nnt
 429       iatel_e=nct-3
 430       do i=iatel_s,iatel_e
 431         ielstart(i)=i+2
 432         ielend(i)=nct-1
 433       enddo
 434 #endif
 435       if (lprint) then
 436         write (iout,'(a)') 'Electrostatic interaction array:'
 437         do i=iatel_s,iatel_e
 438           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
 439         enddo
 440       endif ! lprint
 441 c     iscp=3
 442       iscp=2
 443 C Partition the SC-p interaction array
 444 #ifdef MPL
 445       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
 446       call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
 447       if (lprint)
 448      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 449      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
 450      &               ' my_scp_inde',my_scp_inde
 451       iatscp_s=0
 452       iatscp_e=0
 453       ind_scpint=0
 454       ind_scpint_old=0
 455       do i=nnt,nct-1
 456         if (i.lt.nnt+iscp) then
 457 cd        write (iout,*) 'i.le.nnt+iscp'
 458           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 459      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
 460      &      iscpend(i,1),*14)
 461         else if (i.gt.nct-iscp) then
 462 cd        write (iout,*) 'i.gt.nct-iscp'
 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         else
 467           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 468      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 469      &      iscpend(i,1),*14)
 470           ii=nscp_gr(i)+1
 471           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 472      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
 473      &      iscpend(i,ii),*14)
 474         endif
 475       enddo ! i
 476    14 continue
 477 #else
 478       iatscp_s=nnt
 479       iatscp_e=nct-1
 480       do i=nnt,nct-1
 481         if (i.lt.nnt+iscp) then
 482           nscp_gr(i)=1
 483           iscpstart(i,1)=i+iscp
 484           iscpend(i,1)=nct
 485         elseif (i.gt.nct-iscp) then
 486           nscp_gr(i)=1
 487           iscpstart(i,1)=nnt
 488           iscpend(i,1)=i-iscp
 489         else
 490           nscp_gr(i)=2
 491           iscpstart(i,1)=nnt
 492           iscpend(i,1)=i-iscp
 493           iscpstart(i,2)=i+iscp
 494           iscpend(i,2)=nct
 495         endif
 496       enddo ! i
 497 #endif
 498       if (lprint) then
 499         write (iout,'(a)') 'SC-p interaction array:'
 500         do i=iatscp_s,iatscp_e
 501           write (iout,'(i3,2(2x,2i3))')
 502      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
 503         enddo
 504       endif ! lprint
 505 C Partition local interactions
 506 #ifdef MPL
 507       call int_bounds(nres-2,loc_start,loc_end)
 508       loc_start=loc_start+1
 509       loc_end=loc_end+1
 510       call int_bounds(nres-2,ithet_start,ithet_end)
 511       ithet_start=ithet_start+2
 512       ithet_end=ithet_end+2
 513       call int_bounds(nct-nnt-2,iphi_start,iphi_end)
 514       iphi_start=iphi_start+nnt+2
 515       iphi_end=iphi_end+nnt+2
 516       call int_bounds(nres-3,itau_start,itau_end)
 517       itau_start=itau_start+3
 518       itau_end=itau_end+3
 519       if (lprint) then
 520         write (iout,*) '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         write (*,*) 'Processor:',MyID,
 525      & ' loc_start',loc_start,' loc_end',loc_end,
 526      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 527      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 528       endif
 529       if (fgprocs.gt.1 .and. MyID.eq.BossID) then
 530         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
 531      & nele_int_tot,' electrostatic and ',nscp_int_tot,
 532      & ' SC-p interactions','were distributed among',fgprocs,
 533      & ' fine-grain processors.'
 534       endif
 535 #else
 536       loc_start=2
 537       loc_end=nres-1
 538       ithet_start=3
 539       ithet_end=nres
 540       iphi_start=nnt+3
 541       iphi_end=nct
 542       itau_start=4
 543       itau_end=nres
 544 #endif
 545       return
 546       end
 547 c---------------------------------------------------------------------------
 548       subroutine int_partition(int_index,lower_index,upper_index,atom,
 549      & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
 550       implicit real*8 (a-h,o-z)
 551       include 'DIMENSIONS'
 552       include 'COMMON.IOUNITS'
 553       integer int_index,lower_index,upper_index,atom,at_start,at_end,
 554      & first_atom,last_atom,int_gr,jat_start,jat_end
 555       logical lprn
 556       lprn=.false.
 557       if (lprn) write (iout,*) 'int_index=',int_index
 558       int_index_old=int_index
 559       int_index=int_index+last_atom-first_atom+1
 560       if (lprn)
 561      &   write (iout,*) 'int_index=',int_index,
 562      &               ' int_index_old',int_index_old,
 563      &               ' lower_index=',lower_index,
 564      &               ' upper_index=',upper_index,
 565      &               ' atom=',atom,' first_atom=',first_atom,
 566      &               ' last_atom=',last_atom
 567       if (int_index.ge.lower_index) then
 568         int_gr=int_gr+1
 569         if (at_start.eq.0) then
 570           at_start=atom
 571           jat_start=first_atom-1+lower_index-int_index_old
 572         else
 573           jat_start=first_atom
 574         endif
 575         if (lprn) write (iout,*) 'jat_start',jat_start
 576         if (int_index.ge.upper_index) then
 577           at_end=atom
 578           jat_end=first_atom-1+upper_index-int_index_old
 579           return1
 580         else
 581           jat_end=last_atom
 582         endif
 583         if (lprn) write (iout,*) 'jat_end',jat_end
 584       endif
 585       return
 586       end
 587 c------------------------------------------------------------------------------
 588       subroutine hpb_partition
 589       implicit real*8 (a-h,o-z)
 590       include 'DIMENSIONS'
 591       include 'COMMON.SBRIDGE'
 592       include 'COMMON.IOUNITS'
 593 #ifdef MPL
 594       include 'COMMON.INFO'
 595       call int_bounds(nhpb,link_start,link_end)
 596 #else
 597       link_start=1
 598       link_end=nhpb
 599 #endif
 600 cd    write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 601 cd   &  ' nhpb',nhpb,' link_start=',link_start,
 602 cd   &  ' link_end',link_end
 603       return
 604       end