source/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 'DIMENSIONS.ZSCOPT'
   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.WEIGHTS"
  22       include "COMMON.NAMES"
  23       include "COMMON.TIME1"
  24 C
  25 C The following is just to define auxiliary variables used in angle conversion
  26 C
  27       pi=4.0D0*datan(1.0D0)
  28       dwapi=2.0D0*pi
  29       dwapi3=dwapi/3.0D0
  30       pipol=0.5D0*pi
  31       deg2rad=pi/180.0D0
  32       rad2deg=1.0D0/deg2rad
  33       angmin=10.0D0*deg2rad
  34 C
  35 C Define I/O units.
  36 C
  37       inp=    1
  38       iout=   2
  39       ipdbin= 3
  40       ipdb=   7
  41       imol2=  4
  42       igeom=  8
  43       intin=  9
  44       ithep= 11
  45       irotam=12
  46       itorp= 13
  47       itordp= 23
  48       ielep= 14
  49       isidep=15
  50       isidep1=22
  51       iscpp=25
  52       icbase=16
  53       ifourier=20
  54       istat= 17
  55       ientin=18
  56       ientout=19
  57       ibond=28
  58       isccor=29
  59 C
  60 C WHAM files
  61 C
  62       ihist=30
  63       iweight=31
  64       izsc=32
  65 C
  66 C Set default weights of the energy terms.
  67 C
  68       wlong=1.0D0
  69       welec=1.0D0
  70       wtor =1.0D0
  71       wang =1.0D0
  72       wscloc=1.0D0
  73       wstrain=1.0D0
  74 C
  75 C Zero out tables.
  76 C
  77       ndih_constr=0
  78       do i=1,maxres2
  79         do j=1,3
  80           c(j,i)=0.0D0
  81           dc(j,i)=0.0D0
  82         enddo
  83       enddo
  84       do i=1,maxres
  85         do j=1,3
  86           xloc(j,i)=0.0D0
  87         enddo
  88       enddo
  89       do i=1,ntyp
  90         do j=1,ntyp
  91           aa(i,j)=0.0D0
  92           bb(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           athet(j,i)=0.0D0
 109           bthet(j,i)=0.0D0
 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=1,maxtor
 136         itortyp(i)=0
 137         do j=1,maxtor
 138           do k=1,maxterm
 139             v1(k,j,i)=0.0D0
 140             v2(k,j,i)=0.0D0
 141           enddo
 142         enddo
 143       enddo
 144       do i=1,maxres
 145         itype(i)=0
 146         itel(i)=0
 147       enddo
 148 C Initialize the bridge arrays
 149       ns=0
 150       nss=0
 151       nhpb=0
 152       do i=1,maxss
 153         iss(i)=0
 154       enddo
 155       do i=1,maxdim
 156         dhpb(i)=0.0D0
 157       enddo
 158       do i=1,maxres
 159         ihpb(i)=0
 160         jhpb(i)=0
 161       enddo
 162 C
 163 C Initialize timing.
 164 C
 165       call set_timers
 166 C
 167 C Initialize variables used in minimization.
 168 C
 169 c     maxfun=5000
 170 c     maxit=2000
 171       maxfun=500
 172       maxit=200
 173       tolf=1.0D-2
 174       rtolf=5.0D-4
 175 C
 176 C Initialize the variables responsible for the mode of gradient storage.
 177 C
 178       nfl=0
 179       icg=1
 180       do i=1,14
 181         do j=1,14
 182           if (print_order(i).eq.j) then
 183             iw(print_order(i))=j
 184             goto 1121
 185           endif
 186         enddo
 187 1121    continue
 188       enddo
 189       calc_grad=.false.
 190 C Set timers and counters for the respective routines
 191       t_func = 0.0d0
 192       t_grad = 0.0d0
 193       t_fhel = 0.0d0
 194       t_fbet = 0.0d0
 195       t_ghel = 0.0d0
 196       t_gbet = 0.0d0
 197       t_viol = 0.0d0
 198       t_gviol = 0.0d0
 199       n_func = 0
 200       n_grad = 0
 201       n_fhel = 0
 202       n_fbet = 0
 203       n_ghel = 0
 204       n_gbet = 0
 205       n_viol = 0
 206       n_gviol = 0
 207       n_map = 0
 208 #ifndef SPLITELE
 209       nprint_ene=nprint_ene-1
 210 #endif
 211       return
 212       end
 213 c-------------------------------------------------------------------------
 214       block data nazwy
 215       implicit real*8 (a-h,o-z)
 216       include 'DIMENSIONS'
 217       include 'DIMENSIONS.ZSCOPT'
 218       include 'COMMON.NAMES'
 219       include 'COMMON.WEIGHTS'
 220       include 'COMMON.FFIELD'
 221       data restyp /
 222      &'CYS','MET','PHE','ILE','LEU','VAL','TRP','TYR','ALA','GLY','THR',
 223      &'SER','GLN','ASN','GLU','ASP','HIS','ARG','LYS','PRO','D'/
 224       data onelet /
 225      &'C','M','F','I','L','V','W','Y','A','G','T',
 226      &'S','Q','N','E','D','H','R','K','P','X'/
 227       data potname /'LJ','LJK','BP','GB','GBV'/
 228       data ename /
 229      &   "EVDW SC-SC","EVDW2 SC-p","EES p-p","ECORR4 ","ECORR5 ",
 230      &   "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",
 231      &   "EBE bend","ESC SCloc","ETORS ","ETORSD ","EHPB","EVDWPP",
 232      &   "EVDW2_14","ESTR","ESCCOR","EDIHC","EVDW_T"/
 233       data wname /
 234      &   "WSC","WSCP","WELEC","WCORR","WCORR5","WCORR6","WEL_LOC",
 235      &   "WTURN3","WTURN4","WTURN6","WANG","WSCLOC","WTOR","WTORD",
 236      &   "WHPB","WVDWPP","WSCP14","WBOND","WSCCOR","WDIHC","WSC"/
 237       data ww0 /1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,
 238      &    1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,1.0d0,0.4d0,1.0d0,1.0d0,
 239      &    0.0d0,0.0/
 240       data nprint_ene /21/
 241       data print_order /1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,
 242      &  16,15,17,20,21/
 243       end
 244 c---------------------------------------------------------------------------
 245       subroutine init_int_table
 246       implicit real*8 (a-h,o-z)
 247       include 'DIMENSIONS'
 248       include 'DIMENSIONS.ZSCOPT'
 249 #ifdef MPI
 250       include 'mpif.h'
 251 #endif
 252 #ifdef MP
 253       include 'COMMON.INFO'
 254 #endif
 255       include 'COMMON.CHAIN'
 256       include 'COMMON.INTERACT'
 257       include 'COMMON.LOCAL'
 258       include 'COMMON.SBRIDGE'
 259       include 'COMMON.IOUNITS'
 260       logical scheck,lprint
 261 #ifdef MPL
 262       integer my_sc_int(0:max_fg_Procs-1),my_sc_intt(0:max_fg_Procs),
 263      & my_ele_int(0:max_fg_Procs-1),my_ele_intt(0:max_fg_Procs)
 264 C... Determine the numbers of start and end SC-SC interaction
 265 C... to deal with by current processor.
 266       lprint=.false.
 267       if (lprint)
 268      &write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct
 269       n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss
 270       MyRank=MyID-(MyGroup-1)*fgProcs
 271       call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde)
 272       if (lprint)
 273      &  write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 274      &  ' n_sc_int_tot',n_sc_int_tot,' my_sc_inds=',my_sc_inds,
 275      &  ' my_sc_inde',my_sc_inde
 276       ind_sctint=0
 277       iatsc_s=0
 278       iatsc_e=0
 279 #endif
 280       lprint=.false.
 281       do i=1,maxres
 282         nint_gr(i)=0
 283         nscp_gr(i)=0
 284         do j=1,maxint_gr
 285           istart(i,1)=0
 286           iend(i,1)=0
 287           ielstart(i)=0
 288           ielend(i)=0
 289           iscpstart(i,1)=0
 290           iscpend(i,1)=0
 291         enddo
 292       enddo
 293       ind_scint=0
 294       ind_scint_old=0
 295 cd    write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb',
 296 cd   &   (ihpb(i),jhpb(i),i=1,nss)
 297       do i=nnt,nct-1
 298         scheck=.false.
 299         do ii=1,nss
 300           if (ihpb(ii).eq.i+nres) then
 301             scheck=.true.
 302             jj=jhpb(ii)-nres
 303             goto 10
 304           endif
 305         enddo
 306    10   continue
 307 cd      write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj
 308         if (scheck) then
 309           if (jj.eq.i+1) then
 310 #ifdef MPL
 311             write (iout,*) 'jj=i+1'
 312             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 313      & iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 314 #else
 315             nint_gr(i)=1
 316             istart(i,1)=i+2
 317             iend(i,1)=nct
 318 #endif
 319           else if (jj.eq.nct) then
 320 #ifdef MPL
 321             write (iout,*) 'jj=nct'
 322             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 323      &  iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 324 #else
 325             nint_gr(i)=1
 326             istart(i,1)=i+1
 327             iend(i,1)=nct-1
 328 #endif
 329           else
 330 #ifdef MPL
 331             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 332      & iatsc_s,iatsc_e,i+1,jj-1,nint_gr(i),istart(i,1),iend(i,1),*12)
 333             ii=nint_gr(i)+1
 334             call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 335      & iatsc_s,iatsc_e,jj+1,nct,nint_gr(i),istart(i,ii),iend(i,ii),*12)
 336 #else
 337             nint_gr(i)=2
 338             istart(i,1)=i+1
 339             iend(i,1)=jj-1
 340             istart(i,2)=jj+1
 341             iend(i,2)=nct
 342 #endif
 343           endif
 344         else
 345 #ifdef MPL
 346           call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,
 347      &    iatsc_s,iatsc_e,i+1,nct,nint_gr(i),istart(i,1),iend(i,1),*12)
 348 #else
 349           nint_gr(i)=1
 350           istart(i,1)=i+1
 351           iend(i,1)=nct
 352           ind_scint=int_scint+nct-i
 353 #endif
 354         endif
 355 #ifdef MPL
 356         ind_scint_old=ind_scint
 357 #endif
 358       enddo
 359    12 continue
 360 #ifndef MPL
 361       iatsc_s=nnt
 362       iatsc_e=nct-1
 363 #endif
 364 #ifdef MPL
 365       if (lprint) then
 366         write (iout,*) 'Processor',MyID,' Group',MyGroup
 367         write (iout,*) 'iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e
 368       endif
 369 #endif
 370       if (lprint) then
 371       write (iout,'(a)') 'Interaction array:'
 372       do i=iatsc_s,iatsc_e
 373         write (iout,'(i3,2(2x,2i3))')
 374      & i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i))
 375       enddo
 376       endif
 377       ispp=2
 378 #ifdef MPL
 379 C Now partition the electrostatic-interaction array
 380       npept=nct-nnt
 381       nele_int_tot=(npept-ispp)*(npept-ispp+1)/2
 382       call int_bounds(nele_int_tot,my_ele_inds,my_ele_inde)
 383       if (lprint)
 384      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 385      &  ' nele_int_tot',nele_int_tot,' my_ele_inds=',my_ele_inds,
 386      &               ' my_ele_inde',my_ele_inde
 387       iatel_s=0
 388       iatel_e=0
 389       ind_eleint=0
 390       ind_eleint_old=0
 391       do i=nnt,nct-3
 392         ijunk=0
 393         call int_partition(ind_eleint,my_ele_inds,my_ele_inde,i,
 394      &    iatel_s,iatel_e,i+ispp,nct-1,ijunk,ielstart(i),ielend(i),*13)
 395       enddo ! i
 396    13 continue
 397 #else
 398       iatel_s=nnt
 399       iatel_e=nct-3
 400       do i=iatel_s,iatel_e
 401         ielstart(i)=i+2
 402         ielend(i)=nct-1
 403       enddo
 404 #endif
 405       if (lprint) then
 406         write (iout,'(a)') 'Electrostatic interaction array:'
 407         do i=iatel_s,iatel_e
 408           write (iout,'(i3,2(2x,2i3))') i,ielstart(i),ielend(i)
 409         enddo
 410       endif ! lprint
 411 c     iscp=3
 412       iscp=2
 413 C Partition the SC-p interaction array
 414 #ifdef MPL
 415       nscp_int_tot=(npept-iscp+1)*(npept-iscp+1)
 416       call int_bounds(nscp_int_tot,my_scp_inds,my_scp_inde)
 417       if (lprint)
 418      & write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 419      &  ' nscp_int_tot',nscp_int_tot,' my_scp_inds=',my_scp_inds,
 420      &               ' my_scp_inde',my_scp_inde
 421       iatscp_s=0
 422       iatscp_e=0
 423       ind_scpint=0
 424       ind_scpint_old=0
 425       do i=nnt,nct-1
 426         if (i.lt.nnt+iscp) then
 427 cd        write (iout,*) 'i.le.nnt+iscp'
 428           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 429      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,1),
 430      &      iscpend(i,1),*14)
 431         else if (i.gt.nct-iscp) then
 432 cd        write (iout,*) 'i.gt.nct-iscp'
 433           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 434      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 435      &      iscpend(i,1),*14)
 436         else
 437           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 438      &      iatscp_s,iatscp_e,nnt,i-iscp,nscp_gr(i),iscpstart(i,1),
 439      &      iscpend(i,1),*14)
 440           ii=nscp_gr(i)+1
 441           call int_partition(ind_scpint,my_scp_inds,my_scp_inde,i,
 442      &      iatscp_s,iatscp_e,i+iscp,nct,nscp_gr(i),iscpstart(i,ii),
 443      &      iscpend(i,ii),*14)
 444         endif
 445       enddo ! i
 446    14 continue
 447 #else
 448       iatscp_s=nnt
 449       iatscp_e=nct-1
 450       do i=nnt,nct-1
 451         if (i.lt.nnt+iscp) then
 452           nscp_gr(i)=1
 453           iscpstart(i,1)=i+iscp
 454           iscpend(i,1)=nct
 455         elseif (i.gt.nct-iscp) then
 456           nscp_gr(i)=1
 457           iscpstart(i,1)=nnt
 458           iscpend(i,1)=i-iscp
 459         else
 460           nscp_gr(i)=2
 461           iscpstart(i,1)=nnt
 462           iscpend(i,1)=i-iscp
 463           iscpstart(i,2)=i+iscp
 464           iscpend(i,2)=nct
 465         endif
 466       enddo ! i
 467 #endif
 468       if (lprint) then
 469         write (iout,'(a)') 'SC-p interaction array:'
 470         do i=iatscp_s,iatscp_e
 471           write (iout,'(i3,2(2x,2i3))')
 472      &         i,(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i))
 473         enddo
 474       endif ! lprint
 475 C Partition local interactions
 476 #ifdef MPL
 477       call int_bounds(nres-2,loc_start,loc_end)
 478       loc_start=loc_start+1
 479       loc_end=loc_end+1
 480       call int_bounds(nres-2,ithet_start,ithet_end)
 481       ithet_start=ithet_start+2
 482       ithet_end=ithet_end+2
 483       call int_bounds(nct-nnt-2,iphi_start,iphi_end)
 484       iphi_start=iphi_start+nnt+2
 485       iphi_end=iphi_end+nnt+2
 486       if (lprint) then
 487         write (iout,*) 'Processor:',MyID,
 488      & ' loc_start',loc_start,' loc_end',loc_end,
 489      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 490      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 491         write (*,*) 'Processor:',MyID,
 492      & ' loc_start',loc_start,' loc_end',loc_end,
 493      & ' ithet_start',ithet_start,' ithet_end',ithet_end,
 494      & ' iphi_start',iphi_start,' iphi_end',iphi_end
 495       endif
 496       if (fgprocs.gt.1 .and. MyID.eq.BossID) then
 497         write(iout,'(i10,a,i10,a,i10,a/a,i3,a)') n_sc_int_tot,' SC-SC ',
 498      & nele_int_tot,' electrostatic and ',nscp_int_tot,
 499      & ' SC-p interactions','were distributed among',fgprocs,
 500      & ' fine-grain processors.'
 501       endif
 502 #else
 503       loc_start=2
 504       loc_end=nres-1
 505       ithet_start=3
 506       ithet_end=nres
 507       iphi_start=nnt+3
 508       iphi_end=nct
 509 #endif
 510       return
 511       end
 512 c---------------------------------------------------------------------------
 513       subroutine int_partition(int_index,lower_index,upper_index,atom,
 514      & at_start,at_end,first_atom,last_atom,int_gr,jat_start,jat_end,*)
 515       implicit real*8 (a-h,o-z)
 516       include 'DIMENSIONS'
 517       include 'COMMON.IOUNITS'
 518       integer int_index,lower_index,upper_index,atom,at_start,at_end,
 519      & first_atom,last_atom,int_gr,jat_start,jat_end
 520       logical lprn
 521       lprn=.false.
 522       if (lprn) write (iout,*) 'int_index=',int_index
 523       int_index_old=int_index
 524       int_index=int_index+last_atom-first_atom+1
 525       if (lprn)
 526      &   write (iout,*) 'int_index=',int_index,
 527      &               ' int_index_old',int_index_old,
 528      &               ' lower_index=',lower_index,
 529      &               ' upper_index=',upper_index,
 530      &               ' atom=',atom,' first_atom=',first_atom,
 531      &               ' last_atom=',last_atom
 532       if (int_index.ge.lower_index) then
 533         int_gr=int_gr+1
 534         if (at_start.eq.0) then
 535           at_start=atom
 536           jat_start=first_atom-1+lower_index-int_index_old
 537         else
 538           jat_start=first_atom
 539         endif
 540         if (lprn) write (iout,*) 'jat_start',jat_start
 541         if (int_index.ge.upper_index) then
 542           at_end=atom
 543           jat_end=first_atom-1+upper_index-int_index_old
 544           return1
 545         else
 546           jat_end=last_atom
 547         endif
 548         if (lprn) write (iout,*) 'jat_end',jat_end
 549       endif
 550       return
 551       end
 552 c------------------------------------------------------------------------------
 553       subroutine hpb_partition
 554       implicit real*8 (a-h,o-z)
 555       include 'DIMENSIONS'
 556       include 'COMMON.SBRIDGE'
 557       include 'COMMON.IOUNITS'
 558 #ifdef MPL
 559       include 'COMMON.INFO'
 560       call int_bounds(nhpb,link_start,link_end)
 561 #else
 562       link_start=1
 563       link_end=nhpb
 564 #endif
 565 cd    write (iout,*) 'Processor',MyID,' MyRank',MyRank,
 566 cd   &  ' nhpb',nhpb,' link_start=',link_start,
 567 cd   &  ' link_end',link_end
 568       return
 569       end