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