source/unres/src_MD-M/readrtns_CSA.F

   1       subroutine readrtns
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4 #ifdef MPI
   5       include 'mpif.h'
   6 #endif
   7       include 'COMMON.SETUP'
   8       include 'COMMON.CONTROL'
   9       include 'COMMON.SBRIDGE'
  10       include 'COMMON.IOUNITS'
  11       include 'COMMON.SPLITELE'
  12       logical file_exist
  13 C Read force-field parameters except weights
  14       call parmread
  15 C Read job setup parameters
  16       call read_control
  17 C Read control parameters for energy minimzation if required
  18       if (minim) call read_minim
  19 C Read MCM control parameters if required
  20       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  21 C Read MD control parameters if reqjuired
  22       if (modecalc.eq.12) call read_MDpar
  23 C Read MREMD control parameters if required
  24       if (modecalc.eq.14) then
  25          call read_MDpar
  26          call read_REMDpar
  27       endif
  28 C Read MUCA control parameters if required
  29       if (lmuca) call read_muca
  30 C Read CSA control parameters if required (from fort.40 if exists
  31 C otherwise from general input file)
  32       if (modecalc.eq.8) then
  33        inquire (file="fort.40",exist=file_exist)
  34        if (.not.file_exist) call csaread
  35       endif
  36 cfmc      if (modecalc.eq.10) call mcmfread
  37 C Read molecule information, molecule geometry, energy-term weights, and
  38 C restraints if requested
  39       call molread
  40 C Print restraint information
  41 #ifdef MPI
  42       if (.not. out1file .or. me.eq.king) then
  43 #endif
  44       if (nhpb.gt.nss)
  45      &write (iout,'(a,i5,a)') "The following",nhpb-nss,
  46      & " distance constraints have been imposed"
  47       do i=nss+1,nhpb
  48         write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
  49       enddo
  50 #ifdef MPI
  51       endif
  52 #endif
  53 c      print *,"Processor",myrank," leaves READRTNS"
  54       return
  55       end
  56 C-------------------------------------------------------------------------------
  57       subroutine read_control
  58 C
  59 C Read contorl data
  60 C
  61       implicit real*8 (a-h,o-z)
  62       include 'DIMENSIONS'
  63 #ifdef MP
  64       include 'mpif.h'
  65       logical OKRandom, prng_restart
  66       real*8  r1
  67 #endif
  68       include 'COMMON.IOUNITS'
  69       include 'COMMON.TIME1'
  70       include 'COMMON.THREAD'
  71       include 'COMMON.SBRIDGE'
  72       include 'COMMON.CONTROL'
  73       include 'COMMON.MCM'
  74       include 'COMMON.MAP'
  75       include 'COMMON.HEADER'
  76       include 'COMMON.CSA'
  77       include 'COMMON.CHAIN'
  78       include 'COMMON.MUCA'
  79       include 'COMMON.MD'
  80       include 'COMMON.FFIELD'
  81       include 'COMMON.INTERACT'
  82       include 'COMMON.SETUP'
  83       include 'COMMON.SPLITELE'
  84       include 'COMMON.SHIELD'
  85       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
  86       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
  87       character*80 ucase
  88       character*320 controlcard
  89
  90       nglob_csa=0
  91       eglob_csa=1d99
  92       nmin_csa=0
  93       read (INP,'(a)') titel
  94       call card_concat(controlcard)
  95 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
  96 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
  97       call reada(controlcard,'SEED',seed,0.0D0)
  98       call random_init(seed)
  99 C Set up the time limit (caution! The time must be input in minutes!)
 100       read_cart=index(controlcard,'READ_CART').gt.0
 101       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
 102 C this variable with_theta_constr is the variable which allow to read and execute the
 103 C constrains on theta angles WITH_THETA_CONSTR is the keyword
 104       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
 105       write (iout,*) "with_theta_constr ",with_theta_constr
 106       call readi(controlcard,'SYM',symetr,1)
 107       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
 108       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 109       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 110       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 111       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 112       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 113       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 114       call reada(controlcard,'DRMS',drms,0.1D0)
 115       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 116        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 117        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 118        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 119        write (iout,'(a,f10.1)')'DRMS    = ',drms
 120        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 121        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 122       endif
 123       call readi(controlcard,'NZ_START',nz_start,0)
 124       call readi(controlcard,'NZ_END',nz_end,0)
 125       call readi(controlcard,'IZ_SC',iz_sc,0)
 126       timlim=60.0D0*timlim
 127       safety = 60.0d0*safety
 128       timem=timlim
 129       modecalc=0
 130       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 131       minim=(index(controlcard,'MINIMIZE').gt.0)
 132       dccart=(index(controlcard,'CART').gt.0)
 133       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 134       overlapsc=.not.overlapsc
 135       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
 136       searchsc=.not.searchsc
 137       sideadd=(index(controlcard,'SIDEADD').gt.0)
 138       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 139       outpdb=(index(controlcard,'PDBOUT').gt.0)
 140       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 141       pdbref=(index(controlcard,'PDBREF').gt.0)
 142       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 143       indpdb=index(controlcard,'PDBSTART')
 144       extconf=(index(controlcard,'EXTCONF').gt.0)
 145       AFMlog=(index(controlcard,'AFM'))
 146       selfguide=(index(controlcard,'SELFGUIDE'))
 147       print *,'AFMlog',AFMlog,selfguide,"KUPA"
 148       call readi(controlcard,'IPRINT',iprint,0)
 149 C SHIELD keyword sets if the shielding effect of side-chains is used
 150 C 0 denots no shielding is used all peptide are equally despite the
 151 C solvent accesible area
 152 C 1 the newly introduced function
 153 C 2 reseved for further possible developement
 154       call readi(controlcard,'SHIELD',shield_mode,0)
 155 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 156         write(iout,*) "shield_mode",shield_mode
 157 C      endif
 158       call readi(controlcard,'TORMODE',tor_mode,0)
 159 C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 160         write(iout,*) "torsional and valence angle mode",tor_mode
 161       call readi(controlcard,'MAXGEN',maxgen,10000)
 162       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 163       call readi(controlcard,"KDIAG",kdiag,0)
 164       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
 165       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 166      & write (iout,*) "RESCALE_MODE",rescale_mode
 167       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 168       if (index(controlcard,'REGULAR').gt.0.0D0) then
 169         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 170         modecalc=1
 171         refstr=.true.
 172       endif
 173       if (index(controlcard,'CHECKGRAD').gt.0) then
 174         modecalc=5
 175         if (index(controlcard,'CART').gt.0) then
 176           icheckgrad=1
 177         elseif (index(controlcard,'CARINT').gt.0) then
 178           icheckgrad=2
 179         else
 180           icheckgrad=3
 181         endif
 182       elseif (index(controlcard,'THREAD').gt.0) then
 183         modecalc=2
 184         call readi(controlcard,'THREAD',nthread,0)
 185         if (nthread.gt.0) then
 186           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 187         else
 188           if (fg_rank.eq.0)
 189      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 190           stop 'Error termination in Read_Control.'
 191         endif
 192       else if (index(controlcard,'MCMA').gt.0) then
 193         modecalc=3
 194       else if (index(controlcard,'MCEE').gt.0) then
 195         modecalc=6
 196       else if (index(controlcard,'MULTCONF').gt.0) then
 197         modecalc=4
 198       else if (index(controlcard,'MAP').gt.0) then
 199         modecalc=7
 200         call readi(controlcard,'MAP',nmap,0)
 201       else if (index(controlcard,'CSA').gt.0) then
 202         modecalc=8
 203 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 204 crc
 205 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 206 crc
 207 crc        modecalc=9
 208 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 209 cfmc        modecalc=10
 210       else if (index(controlcard,'SOFTREG').gt.0) then
 211         modecalc=11
 212       else if (index(controlcard,'CHECK_BOND').gt.0) then
 213         modecalc=-1
 214       else if (index(controlcard,'TEST').gt.0) then
 215         modecalc=-2
 216       else if (index(controlcard,'MD').gt.0) then
 217         modecalc=12
 218       else if (index(controlcard,'RE ').gt.0) then
 219         modecalc=14
 220       endif
 221
 222       lmuca=index(controlcard,'MUCA').gt.0
 223       call readi(controlcard,'MUCADYN',mucadyn,0)
 224       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 225       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 226      & then
 227        write (iout,*) 'MUCADYN=',mucadyn
 228        write (iout,*) 'MUCASMOOTH=',muca_smooth
 229       endif
 230
 231       iscode=index(controlcard,'ONE_LETTER')
 232       indphi=index(controlcard,'PHI')
 233       indback=index(controlcard,'BACK')
 234       iranconf=index(controlcard,'RAND_CONF')
 235       i2ndstr=index(controlcard,'USE_SEC_PRED')
 236       gradout=index(controlcard,'GRADOUT').gt.0
 237       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 238 C  DISTCHAINMAX become obsolete for periodic boundry condition
 239       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
 240 C Reading the dimensions of box in x,y,z coordinates
 241       call reada(controlcard,'BOXX',boxxsize,100.0d0)
 242       call reada(controlcard,'BOXY',boxysize,100.0d0)
 243       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
 244 c Cutoff range for interactions
 245       call reada(controlcard,"R_CUT",r_cut,15.0d0)
 246       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 247       call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
 248       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
 249       if (lipthick.gt.0.0d0) then
 250        bordliptop=(boxzsize+lipthick)/2.0
 251        bordlipbot=bordliptop-lipthick
 252 C      endif
 253       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
 254      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
 255       buflipbot=bordlipbot+lipbufthick
 256       bufliptop=bordliptop-lipbufthick
 257       if ((lipbufthick*2.0d0).gt.lipthick)
 258      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
 259       endif
 260       write(iout,*) "bordliptop=",bordliptop
 261       write(iout,*) "bordlipbot=",bordlipbot
 262       write(iout,*) "bufliptop=",bufliptop
 263       write(iout,*) "buflipbot=",buflipbot
 264       write (iout,*) "SHIELD MODE",shield_mode
 265       if (shield_mode.gt.0) then
 266       pi=3.141592d0
 267 C VSolvSphere the volume of solving sphere
 268 C      print *,pi,"pi"
 269 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
 270 C there will be no distinction between proline peptide group and normal peptide
 271 C group in case of shielding parameters
 272       VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
 273       VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
 274       write (iout,*) VSolvSphere,VSolvSphere_div
 275 C long axis of side chain
 276       do i=1,ntyp
 277       long_r_sidechain(i)=vbldsc0(1,i)
 278       short_r_sidechain(i)=sigma0(i)
 279       enddo
 280       buff_shield=1.0d0
 281       endif
 282       if (me.eq.king .or. .not.out1file )
 283      &  write (iout,*) "DISTCHAINMAX",distchainmax
 284
 285       if(me.eq.king.or..not.out1file)
 286      & write (iout,'(2a)') diagmeth(kdiag),
 287      &  ' routine used to diagonalize matrices.'
 288       return
 289       end
 290 c--------------------------------------------------------------------------
 291       subroutine read_REMDpar
 292 C
 293 C Read REMD settings
 294 C
 295       implicit real*8 (a-h,o-z)
 296       include 'DIMENSIONS'
 297       include 'COMMON.IOUNITS'
 298       include 'COMMON.TIME1'
 299       include 'COMMON.MD'
 300 #ifndef LANG0
 301       include 'COMMON.LANGEVIN'
 302 #else
 303       include 'COMMON.LANGEVIN.lang0'
 304 #endif
 305       include 'COMMON.INTERACT'
 306       include 'COMMON.NAMES'
 307       include 'COMMON.GEO'
 308       include 'COMMON.REMD'
 309       include 'COMMON.CONTROL'
 310       include 'COMMON.SETUP'
 311       character*80 ucase
 312       character*320 controlcard
 313       character*3200 controlcard1
 314       integer iremd_m_total
 315
 316       if(me.eq.king.or..not.out1file)
 317      & write (iout,*) "REMD setup"
 318
 319       call card_concat(controlcard)
 320       call readi(controlcard,"NREP",nrep,3)
 321       call readi(controlcard,"NSTEX",nstex,1000)
 322       call reada(controlcard,"RETMIN",retmin,10.0d0)
 323       call reada(controlcard,"RETMAX",retmax,1000.0d0)
 324       mremdsync=(index(controlcard,'SYNC').gt.0)
 325       call readi(controlcard,"NSYN",i_sync_step,100)
 326       restart1file=(index(controlcard,'REST1FILE').gt.0)
 327       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
 328       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
 329       if(max_cache_traj_use.gt.max_cache_traj)
 330      &           max_cache_traj_use=max_cache_traj
 331       if(me.eq.king.or..not.out1file) then
 332 cd       if (traj1file) then
 333 crc caching is in testing - NTWX is not ignored
 334 cd        write (iout,*) "NTWX value is ignored"
 335 cd        write (iout,*) "  trajectory is stored to one file by master"
 336 cd        write (iout,*) "  before exchange at NSTEX intervals"
 337 cd       endif
 338        write (iout,*) "NREP= ",nrep
 339        write (iout,*) "NSTEX= ",nstex
 340        write (iout,*) "SYNC= ",mremdsync
 341        write (iout,*) "NSYN= ",i_sync_step
 342        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
 343       endif
 344       remd_tlist=.false.
 345       if (index(controlcard,'TLIST').gt.0) then
 346          remd_tlist=.true.
 347          call card_concat(controlcard1)
 348          read(controlcard1,*) (remd_t(i),i=1,nrep)
 349          if(me.eq.king.or..not.out1file)
 350      &    write (iout,*)'tlist',(remd_t(i),i=1,nrep)
 351       endif
 352       remd_mlist=.false.
 353       if (index(controlcard,'MLIST').gt.0) then
 354          remd_mlist=.true.
 355          call card_concat(controlcard1)
 356          read(controlcard1,*) (remd_m(i),i=1,nrep)
 357          if(me.eq.king.or..not.out1file) then
 358           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
 359           iremd_m_total=0
 360           do i=1,nrep
 361            iremd_m_total=iremd_m_total+remd_m(i)
 362           enddo
 363            write (iout,*) 'Total number of replicas ',iremd_m_total
 364           endif
 365          endif
 366       if(me.eq.king.or..not.out1file)
 367      &   write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
 368       return
 369       end
 370 c--------------------------------------------------------------------------
 371       subroutine read_MDpar
 372 C
 373 C Read MD settings
 374 C
 375       implicit real*8 (a-h,o-z)
 376       include 'DIMENSIONS'
 377       include 'COMMON.IOUNITS'
 378       include 'COMMON.TIME1'
 379       include 'COMMON.MD'
 380 #ifndef LANG0
 381       include 'COMMON.LANGEVIN'
 382 #else
 383       include 'COMMON.LANGEVIN.lang0'
 384 #endif
 385       include 'COMMON.INTERACT'
 386       include 'COMMON.NAMES'
 387       include 'COMMON.GEO'
 388       include 'COMMON.SETUP'
 389       include 'COMMON.CONTROL'
 390       include 'COMMON.SPLITELE'
 391       character*80 ucase
 392       character*320 controlcard
 393
 394       call card_concat(controlcard)
 395       call readi(controlcard,"NSTEP",n_timestep,1000000)
 396       call readi(controlcard,"NTWE",ntwe,100)
 397       call readi(controlcard,"NTWX",ntwx,1000)
 398       call reada(controlcard,"DT",d_time,1.0d-1)
 399       call reada(controlcard,"DVMAX",dvmax,2.0d1)
 400       call reada(controlcard,"DAMAX",damax,1.0d1)
 401       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
 402       call readi(controlcard,"LANG",lang,0)
 403       RESPA = index(controlcard,"RESPA") .gt. 0
 404       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
 405       ntime_split0=ntime_split
 406       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
 407       ntime_split0=ntime_split
 408 c      call reada(controlcard,"R_CUT",r_cut,2.0d0)
 409 c      call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 410       rest = index(controlcard,"REST").gt.0
 411       tbf = index(controlcard,"TBF").gt.0
 412       usampl = index(controlcard,"USAMPL").gt.0
 413
 414       mdpdb = index(controlcard,"MDPDB").gt.0
 415       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 416       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
 417       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
 418       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
 419       if (count_reset_moment.eq.0) count_reset_moment=1000000000
 420       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
 421       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
 422       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
 423       if (count_reset_vel.eq.0) count_reset_vel=1000000000
 424       large = index(controlcard,"LARGE").gt.0
 425       print_compon = index(controlcard,"PRINT_COMPON").gt.0
 426       rattle = index(controlcard,"RATTLE").gt.0
 427 c  if performing umbrella sampling, fragments constrained are read from the fragment file
 428       nset=0
 429       if(usampl) then
 430         call read_fragments
 431       endif
 432
 433       if(me.eq.king.or..not.out1file) then
 434        write (iout,*)
 435        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
 436        write (iout,*)
 437        write (iout,'(a)') "The units are:"
 438        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
 439        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
 440      &  " acceleration: angstrom/(48.9 fs)**2"
 441        write (iout,'(a)') "energy: kcal/mol, temperature: K"
 442        write (iout,*)
 443        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
 444        write (iout,'(a60,f10.5,a)')
 445      &  "Initial time step of numerical integration:",d_time,
 446      &  " natural units"
 447        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
 448        if (RESPA) then
 449         write (iout,'(2a,i4,a)')
 450      &    "A-MTS algorithm used; initial time step for fast-varying",
 451      &    " short-range forces split into",ntime_split," steps."
 452         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
 453      &   r_cut," lambda",rlamb
 454        endif
 455        write (iout,'(2a,f10.5)')
 456      &  "Maximum acceleration threshold to reduce the time step",
 457      &  "/increase split number:",damax
 458        write (iout,'(2a,f10.5)')
 459      &  "Maximum predicted energy drift to reduce the timestep",
 460      &  "/increase split number:",edriftmax
 461        write (iout,'(a60,f10.5)')
 462      & "Maximum velocity threshold to reduce velocities:",dvmax
 463        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
 464        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
 465        if (rattle) write (iout,'(a60)')
 466      &  "Rattle algorithm used to constrain the virtual bonds"
 467       endif
 468       reset_fricmat=1000
 469       if (lang.gt.0) then
 470         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
 471         call reada(controlcard,"RWAT",rwat,1.4d0)
 472         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
 473         surfarea=index(controlcard,"SURFAREA").gt.0
 474         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
 475         if(me.eq.king.or..not.out1file)then
 476          write (iout,'(/a,$)') "Langevin dynamics calculation"
 477          if (lang.eq.1) then
 478           write (iout,'(a/)')
 479      &      " with direct integration of Langevin equations"
 480          else if (lang.eq.2) then
 481           write (iout,'(a/)') " with TINKER stochasic MD integrator"
 482          else if (lang.eq.3) then
 483           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
 484          else if (lang.eq.4) then
 485           write (iout,'(a/)') " in overdamped mode"
 486          else
 487           write (iout,'(//a,i5)')
 488      &      "=========== ERROR: Unknown Langevin dynamics mode:",lang
 489           stop
 490          endif
 491          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 492          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
 493          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
 494          write (iout,'(a60,f10.5)')
 495      &   "Scaling factor of the friction forces:",scal_fric
 496          if (surfarea) write (iout,'(2a,i10,a)')
 497      &     "Friction coefficients will be scaled by solvent-accessible",
 498      &     " surface area every",reset_fricmat," steps."
 499         endif
 500 c Calculate friction coefficients and bounds of stochastic forces
 501         eta=6*pi*cPoise*etawat
 502         if(me.eq.king.or..not.out1file)
 503      &   write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
 504      &   ,eta
 505         gamp=scal_fric*(pstok+rwat)*eta
 506         stdfp=dsqrt(2*Rb*t_bath/d_time)
 507         do i=1,ntyp
 508           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
 509           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 510         enddo
 511         if(me.eq.king.or..not.out1file)then
 512          write (iout,'(/2a/)')
 513      &   "Radii of site types and friction coefficients and std's of",
 514      &   " stochastic forces of fully exposed sites"
 515          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
 516          do i=1,ntyp
 517           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
 518      &     gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
 519          enddo
 520         endif
 521       else if (tbf) then
 522         if(me.eq.king.or..not.out1file)then
 523          write (iout,'(a)') "Berendsen bath calculation"
 524          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 525          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
 526          if (reset_moment)
 527      &   write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
 528      &   count_reset_moment," steps"
 529          if (reset_vel)
 530      &    write (iout,'(a,i10,a)')
 531      &    "Velocities will be reset at random every",count_reset_vel,
 532      &   " steps"
 533         endif
 534       else
 535         if(me.eq.king.or..not.out1file)
 536      &   write (iout,'(a31)') "Microcanonical mode calculation"
 537       endif
 538       if(me.eq.king.or..not.out1file)then
 539        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
 540        if (usampl) then
 541           write(iout,*) "MD running with constraints."
 542           write(iout,*) "Equilibration time ", eq_time, " mtus."
 543           write(iout,*) "Constraining ", nfrag," fragments."
 544           write(iout,*) "Length of each fragment, weight and q0:"
 545           do iset=1,nset
 546            write (iout,*) "Set of restraints #",iset
 547            do i=1,nfrag
 548               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
 549      &           ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
 550            enddo
 551            write(iout,*) "constraints between ", npair, "fragments."
 552            write(iout,*) "constraint pairs, weights and q0:"
 553            do i=1,npair
 554             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
 555      &             ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
 556            enddo
 557            write(iout,*) "angle constraints within ", nfrag_back,
 558      &      "backbone fragments."
 559            write(iout,*) "fragment, weights:"
 560            do i=1,nfrag_back
 561             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
 562      &         ifrag_back(2,i,iset),wfrag_back(1,i,iset),
 563      &         wfrag_back(2,i,iset),wfrag_back(3,i,iset)
 564            enddo
 565           enddo
 566         iset=mod(kolor,nset)+1
 567        endif
 568       endif
 569       if(me.eq.king.or..not.out1file)
 570      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
 571       return
 572       end
 573 c------------------------------------------------------------------------------
 574       subroutine molread
 575 C
 576 C Read molecular data.
 577 C
 578       implicit real*8 (a-h,o-z)
 579       include 'DIMENSIONS'
 580 #ifdef MPI
 581       include 'mpif.h'
 582       integer error_msg
 583 #endif
 584       include 'COMMON.IOUNITS'
 585       include 'COMMON.GEO'
 586       include 'COMMON.VAR'
 587       include 'COMMON.INTERACT'
 588       include 'COMMON.LOCAL'
 589       include 'COMMON.NAMES'
 590       include 'COMMON.CHAIN'
 591       include 'COMMON.FFIELD'
 592       include 'COMMON.SBRIDGE'
 593       include 'COMMON.HEADER'
 594       include 'COMMON.CONTROL'
 595       include 'COMMON.DBASE'
 596       include 'COMMON.THREAD'
 597       include 'COMMON.CONTACTS'
 598       include 'COMMON.TORCNSTR'
 599       include 'COMMON.TIME1'
 600       include 'COMMON.BOUNDS'
 601       include 'COMMON.MD'
 602       include 'COMMON.SETUP'
 603       include 'COMMON.SHIELD'
 604       character*4 sequence(maxres)
 605       integer rescode
 606       double precision x(maxvar)
 607       character*256 pdbfile
 608       character*400 weightcard
 609       character*80 weightcard_t,ucase
 610       dimension itype_pdb(maxres)
 611       common /pizda/ itype_pdb
 612       logical seq_comp,fail
 613       double precision energia(0:n_ene)
 614       integer ilen
 615       external ilen
 616 C
 617 C Body
 618 C
 619 C Read weights of the subsequent energy terms.
 620        call card_concat(weightcard)
 621        call reada(weightcard,'WLONG',wlong,1.0D0)
 622        call reada(weightcard,'WSC',wsc,wlong)
 623        call reada(weightcard,'WSCP',wscp,wlong)
 624        call reada(weightcard,'WELEC',welec,1.0D0)
 625        call reada(weightcard,'WVDWPP',wvdwpp,welec)
 626        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 627        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 628        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 629        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 630        call reada(weightcard,'WTURN3',wturn3,1.0D0)
 631        call reada(weightcard,'WTURN4',wturn4,1.0D0)
 632        call reada(weightcard,'WTURN6',wturn6,1.0D0)
 633        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 634        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 635        call reada(weightcard,'WBOND',wbond,1.0D0)
 636        call reada(weightcard,'WTOR',wtor,1.0D0)
 637        call reada(weightcard,'WTORD',wtor_d,1.0D0)
 638        call reada(weightcard,'WANG',wang,1.0D0)
 639        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 640        call reada(weightcard,'SCAL14',scal14,0.4D0)
 641        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 642        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 643        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 644        call reada(weightcard,'TEMP0',temp0,300.0d0)
 645        call reada(weightcard,'WSHIELD',wshield,1.0d0)
 646        call reada(weightcard,'WLT',wliptran,0.0D0)
 647        if (index(weightcard,'SOFT').gt.0) ipot=6
 648 C 12/1/95 Added weight for the multi-body term WCORR
 649        call reada(weightcard,'WCORRH',wcorr,1.0D0)
 650        if (wcorr4.gt.0.0d0) wcorr=wcorr4
 651        weights(1)=wsc
 652        weights(2)=wscp
 653        weights(3)=welec
 654        weights(4)=wcorr
 655        weights(5)=wcorr5
 656        weights(6)=wcorr6
 657        weights(7)=wel_loc
 658        weights(8)=wturn3
 659        weights(9)=wturn4
 660        weights(10)=wturn6
 661        weights(11)=wang
 662        weights(12)=wscloc
 663        weights(13)=wtor
 664        weights(14)=wtor_d
 665        weights(15)=wstrain
 666        weights(16)=wvdwpp
 667        weights(17)=wbond
 668        weights(18)=scal14
 669        weights(21)=wsccor
 670       if(me.eq.king.or..not.out1file)
 671      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 672      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 673      &  wturn4,wturn6
 674    10 format (/'Energy-term weights (unscaled):'//
 675      & 'WSCC=   ',f10.6,' (SC-SC)'/
 676      & 'WSCP=   ',f10.6,' (SC-p)'/
 677      & 'WELEC=  ',f10.6,' (p-p electr)'/
 678      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 679      & 'WBOND=  ',f10.6,' (stretching)'/
 680      & 'WANG=   ',f10.6,' (bending)'/
 681      & 'WSCLOC= ',f10.6,' (SC local)'/
 682      & 'WTOR=   ',f10.6,' (torsional)'/
 683      & 'WTORD=  ',f10.6,' (double torsional)'/
 684      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 685      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 686      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 687      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 688      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 689      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
 690      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 691      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 692      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 693       if(me.eq.king.or..not.out1file)then
 694        if (wcorr4.gt.0.0d0) then
 695         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 696      &   'between contact pairs of peptide groups'
 697         write (iout,'(2(a,f5.3/))')
 698      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 699      &  'Range of quenching the correlation terms:',2*delt_corr
 700        else if (wcorr.gt.0.0d0) then
 701         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 702      &   'between contact pairs of peptide groups'
 703        endif
 704        write (iout,'(a,f8.3)')
 705      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 706        write (iout,'(a,f8.3)')
 707      &  'General scaling factor of SC-p interactions:',scalscp
 708       endif
 709       r0_corr=cutoff_corr-delt_corr
 710       do i=1,ntyp
 711         aad(i,1)=scalscp*aad(i,1)
 712         aad(i,2)=scalscp*aad(i,2)
 713         bad(i,1)=scalscp*bad(i,1)
 714         bad(i,2)=scalscp*bad(i,2)
 715       enddo
 716       call rescale_weights(t_bath)
 717       if(me.eq.king.or..not.out1file)
 718      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 719      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 720      &  wturn4,wturn6
 721    22 format (/'Energy-term weights (scaled):'//
 722      & 'WSCC=   ',f10.6,' (SC-SC)'/
 723      & 'WSCP=   ',f10.6,' (SC-p)'/
 724      & 'WELEC=  ',f10.6,' (p-p electr)'/
 725      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 726      & 'WBOND=  ',f10.6,' (stretching)'/
 727      & 'WANG=   ',f10.6,' (bending)'/
 728      & 'WSCLOC= ',f10.6,' (SC local)'/
 729      & 'WTOR=   ',f10.6,' (torsional)'/
 730      & 'WTORD=  ',f10.6,' (double torsional)'/
 731      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 732      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 733      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 734      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 735      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 736      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
 737      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 738      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 739      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 740       if(me.eq.king.or..not.out1file)
 741      & write (iout,*) "Reference temperature for weights calculation:",
 742      &  temp0
 743       call reada(weightcard,"D0CM",d0cm,3.78d0)
 744       call reada(weightcard,"AKCM",akcm,15.1d0)
 745       call reada(weightcard,"AKTH",akth,11.0d0)
 746       call reada(weightcard,"AKCT",akct,12.0d0)
 747       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 748       call reada(weightcard,"V2SS",v2ss,7.61d0)
 749       call reada(weightcard,"V3SS",v3ss,13.7d0)
 750       call reada(weightcard,"EBR",ebr,-5.50D0)
 751       call reada(weightcard,"ATRISS",atriss,0.301D0)
 752       call reada(weightcard,"BTRISS",btriss,0.021D0)
 753       call reada(weightcard,"CTRISS",ctriss,1.001D0)
 754       call reada(weightcard,"DTRISS",dtriss,1.001D0)
 755       write (iout,*) "ATRISS=", atriss
 756       write (iout,*) "BTRISS=", btriss
 757       write (iout,*) "CTRISS=", ctriss
 758       write (iout,*) "DTRISS=", dtriss
 759       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 760       do i=1,maxres
 761         dyn_ss_mask(i)=.false.
 762       enddo
 763       do i=1,maxres-1
 764         do j=i+1,maxres
 765           dyn_ssbond_ij(i,j)=1.0d300
 766         enddo
 767       enddo
 768       call reada(weightcard,"HT",Ht,0.0D0)
 769       if (dyn_ss) then
 770         ss_depth=ebr/wsc-0.25*eps(1,1)
 771         Ht=Ht/wsc-0.25*eps(1,1)
 772         akcm=akcm*wstrain/wsc
 773         akth=akth*wstrain/wsc
 774         akct=akct*wstrain/wsc
 775         v1ss=v1ss*wstrain/wsc
 776         v2ss=v2ss*wstrain/wsc
 777         v3ss=v3ss*wstrain/wsc
 778       else
 779         if (wstrain.ne.0.0) then
 780          ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 781         else
 782           ss_depth=0.0
 783         endif
 784       endif
 785
 786       if(me.eq.king.or..not.out1file) then
 787        write (iout,*) "Parameters of the SS-bond potential:"
 788        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 789      & " AKCT",akct
 790        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 791        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 792        write (iout,*)" HT",Ht
 793        print *,'indpdb=',indpdb,' pdbref=',pdbref
 794       endif
 795       if (indpdb.gt.0 .or. pdbref) then
 796         read(inp,'(a)') pdbfile
 797         if(me.eq.king.or..not.out1file)
 798      &   write (iout,'(2a)') 'PDB data will be read from file ',
 799      &   pdbfile(:ilen(pdbfile))
 800         open(ipdbin,file=pdbfile,status='old',err=33)
 801         goto 34
 802   33    write (iout,'(a)') 'Error opening PDB file.'
 803         stop
 804   34    continue
 805 c        write (iout,*) 'Begin reading pdb data'
 806 c        call flush(iout)
 807         call readpdb
 808 c        write (iout,*) 'Finished reading pdb data'
 809 c        call flush(iout)
 810         if(me.eq.king.or..not.out1file)
 811      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 812      &   ' nstart_sup=',nstart_sup
 813         do i=1,nres
 814           itype_pdb(i)=itype(i)
 815         enddo
 816         close (ipdbin)
 817         nnt=nstart_sup
 818         nct=nstart_sup+nsup-1
 819         call contact(.false.,ncont_ref,icont_ref,co)
 820
 821         if (sideadd) then
 822          if(me.eq.king.or..not.out1file)
 823      &    write(iout,*)'Adding sidechains'
 824          maxsi=1000
 825          do i=2,nres-1
 826           iti=itype(i)
 827           if (iti.ne.10 .and. itype(i).ne.ntyp1) then
 828             nsi=0
 829             fail=.true.
 830             do while (fail.and.nsi.le.maxsi)
 831               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 832               nsi=nsi+1
 833             enddo
 834             if(fail) write(iout,*)'Adding sidechain failed for res ',
 835      &              i,' after ',nsi,' trials'
 836           endif
 837          enddo
 838         endif
 839       endif
 840       if (indpdb.eq.0) then
 841 C Read sequence if not taken from the pdb file.
 842         read (inp,*) nres
 843 c        print *,'nres=',nres
 844         if (iscode.gt.0) then
 845           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 846         else
 847           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 848         endif
 849 C Convert sequence to numeric code
 850         do i=1,nres
 851           itype(i)=rescode(i,sequence(i),iscode)
 852         enddo
 853 C Assign initial virtual bond lengths
 854         do i=2,nres
 855           vbld(i)=vbl
 856           vbld_inv(i)=vblinv
 857         enddo
 858         do i=2,nres-1
 859           vbld(i+nres)=dsc(iabs(itype(i)))
 860           vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
 861 c          write (iout,*) "i",i," itype",itype(i),
 862 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 863         enddo
 864       endif
 865 c      print *,nres
 866 c      print '(20i4)',(itype(i),i=1,nres)
 867       do i=1,nres
 868 #ifdef PROCOR
 869         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 870 #else
 871         if (itype(i).eq.ntyp1) then
 872 #endif
 873           itel(i)=0
 874 #ifdef PROCOR
 875         else if (iabs(itype(i+1)).ne.20) then
 876 #else
 877         else if (iabs(itype(i)).ne.20) then
 878 #endif
 879           itel(i)=1
 880         else
 881           itel(i)=2
 882         endif
 883       enddo
 884       if(me.eq.king.or..not.out1file)then
 885        write (iout,*) "ITEL"
 886        do i=1,nres-1
 887          write (iout,*) i,itype(i),itel(i)
 888        enddo
 889        print *,'Call Read_Bridge.'
 890       endif
 891       call read_bridge
 892 C 8/13/98 Set limits to generating the dihedral angles
 893       do i=1,nres
 894         phibound(1,i)=-pi
 895         phibound(2,i)=pi
 896       enddo
 897       read (inp,*) ndih_constr
 898       if (ndih_constr.gt.0) then
 899 C        read (inp,*) ftors
 900         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
 901      &  i=1,ndih_constr)
 902         if(me.eq.king.or..not.out1file)then
 903          write (iout,*)
 904      &   'There are',ndih_constr,' constraints on phi angles.'
 905          do i=1,ndih_constr
 906           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
 907      &    ftors(i)
 908          enddo
 909         endif
 910         do i=1,ndih_constr
 911           phi0(i)=deg2rad*phi0(i)
 912           drange(i)=deg2rad*drange(i)
 913         enddo
 914 C        if(me.eq.king.or..not.out1file)
 915 C     &   write (iout,*) 'FTORS',ftors
 916         do i=1,ndih_constr
 917           ii = idih_constr(i)
 918           phibound(1,ii) = phi0(i)-drange(i)
 919           phibound(2,ii) = phi0(i)+drange(i)
 920         enddo
 921       endif
 922 C first setting the theta boundaries to 0 to pi
 923 C this mean that there is no energy penalty for any angle occuring this can be applied
 924 C for generate random conformation but is not implemented in this way
 925 C      do i=1,nres
 926 C        thetabound(1,i)=0
 927 C        thetabound(2,i)=pi
 928 C      enddo
 929 C begin reading theta constrains this is quartic constrains allowing to
 930 C have smooth second derivative
 931       if (with_theta_constr) then
 932 C with_theta_constr is keyword allowing for occurance of theta constrains
 933       read (inp,*) ntheta_constr
 934 C ntheta_constr is the number of theta constrains
 935       if (ntheta_constr.gt.0) then
 936 C        read (inp,*) ftors
 937         read (inp,*) (itheta_constr(i),theta_constr0(i),
 938      &  theta_drange(i),for_thet_constr(i),
 939      &  i=1,ntheta_constr)
 940 C the above code reads from 1 to ntheta_constr
 941 C itheta_constr(i) residue i for which is theta_constr
 942 C theta_constr0 the global minimum value
 943 C theta_drange is range for which there is no energy penalty
 944 C for_thet_constr is the force constant for quartic energy penalty
 945 C E=k*x**4
 946         if(me.eq.king.or..not.out1file)then
 947          write (iout,*)
 948      &   'There are',ntheta_constr,' constraints on phi angles.'
 949          do i=1,ntheta_constr
 950           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
 951      &    theta_drange(i),
 952      &    for_thet_constr(i)
 953          enddo
 954         endif
 955         do i=1,ntheta_constr
 956           theta_constr0(i)=deg2rad*theta_constr0(i)
 957           theta_drange(i)=deg2rad*theta_drange(i)
 958         enddo
 959 C        if(me.eq.king.or..not.out1file)
 960 C     &   write (iout,*) 'FTORS',ftors
 961 C        do i=1,ntheta_constr
 962 C          ii = itheta_constr(i)
 963 C          thetabound(1,ii) = phi0(i)-drange(i)
 964 C          thetabound(2,ii) = phi0(i)+drange(i)
 965 C        enddo
 966       endif ! ntheta_constr.gt.0
 967       endif! with_theta_constr
 968 C
 969 C      with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
 970 C      write (iout,*) "with_dihed_constr ",with_dihed_constr
 971       nnt=1
 972 #ifdef MPI
 973       if (me.eq.king) then
 974 #endif
 975        write (iout,'(a)') 'Boundaries in phi angle sampling:'
 976        do i=1,nres
 977          write (iout,'(a3,i5,2f10.1)')
 978      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 979        enddo
 980 #ifdef MP
 981       endif
 982 #endif
 983       nct=nres
 984 cd      print *,'NNT=',NNT,' NCT=',NCT
 985       if (itype(1).eq.ntyp1) nnt=2
 986       if (itype(nres).eq.ntyp1) nct=nct-1
 987       if (pdbref) then
 988         if(me.eq.king.or..not.out1file)
 989      &   write (iout,'(a,i3)') 'nsup=',nsup
 990         nstart_seq=nnt
 991         if (nsup.le.(nct-nnt+1)) then
 992           do i=0,nct-nnt+1-nsup
 993             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 994               nstart_seq=nnt+i
 995               goto 111
 996             endif
 997           enddo
 998           write (iout,'(a)')
 999      &            'Error - sequences to be superposed do not match.'
1000           stop
1001         else
1002           do i=0,nsup-(nct-nnt+1)
1003             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1004      &      then
1005               nstart_sup=nstart_sup+i
1006               nsup=nct-nnt+1
1007               goto 111
1008             endif
1009           enddo
1010           write (iout,'(a)')
1011      &            'Error - sequences to be superposed do not match.'
1012         endif
1013   111   continue
1014         if (nsup.eq.0) nsup=nct-nnt
1015         if (nstart_sup.eq.0) nstart_sup=nnt
1016         if (nstart_seq.eq.0) nstart_seq=nnt
1017         if(me.eq.king.or..not.out1file)
1018      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1019      &                 ' nstart_seq=',nstart_seq
1020       endif
1021 c--- Zscore rms -------
1022       if (nz_start.eq.0) nz_start=nnt
1023       if (nz_end.eq.0 .and. nsup.gt.0) then
1024         nz_end=nnt+nsup-1
1025       else if (nz_end.eq.0) then
1026         nz_end=nct
1027       endif
1028       if(me.eq.king.or..not.out1file)then
1029        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1030        write (iout,*) 'IZ_SC=',iz_sc
1031       endif
1032 c----------------------
1033       call init_int_table
1034       if (refstr) then
1035         if (.not.pdbref) then
1036           call read_angles(inp,*38)
1037           goto 39
1038    38     write (iout,'(a)') 'Error reading reference structure.'
1039 #ifdef MPI
1040           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1041           stop 'Error reading reference structure'
1042 #endif
1043    39     call chainbuild
1044           call setup_var
1045 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1046           nstart_sup=nnt
1047           nstart_seq=nnt
1048           nsup=nct-nnt+1
1049           kkk=1
1050           do i=1,2*nres
1051             do j=1,3
1052               cref(j,i,kkk)=c(j,i)
1053             enddo
1054           enddo
1055           call contact(.true.,ncont_ref,icont_ref,co)
1056         endif
1057         endif
1058         print *, "A TU"
1059         write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1060         call flush(iout)
1061         if (constr_dist.gt.0) call read_dist_constr
1062         write (iout,*) "After read_dist_constr nhpb",nhpb
1063         if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1064         call hpb_partition
1065         if(me.eq.king.or..not.out1file)
1066      &   write (iout,*) 'Contact order:',co
1067         if (pdbref) then
1068         if(me.eq.king.or..not.out1file)
1069      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1070         do i=1,ncont_ref
1071           do j=1,2
1072             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1073           enddo
1074           if(me.eq.king.or..not.out1file)
1075      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1076      &     icont_ref(1,i),' ',
1077      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1078         enddo
1079         endif
1080 C      endif
1081       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1082      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1083      &    modecalc.ne.10) then
1084 C If input structure hasn't been supplied from the PDB file read or generate
1085 C initial geometry.
1086         if (iranconf.eq.0 .and. .not. extconf) then
1087           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1088      &     write (iout,'(a)') 'Initial geometry will be read in.'
1089           if (read_cart) then
1090             read(inp,'(8f10.5)',end=36,err=36)
1091      &       ((c(l,k),l=1,3),k=1,nres),
1092      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1093             write (iout,*) "Exit READ_CART"
1094             write (iout,'(8f10.5)')
1095      &       ((c(l,k),l=1,3),k=1,nres),
1096      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1097             call int_from_cart1(.true.)
1098             write (iout,*) "Finish INT_TO_CART"
1099             do i=1,nres-1
1100               do j=1,3
1101                 dc(j,i)=c(j,i+1)-c(j,i)
1102                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1103               enddo
1104             enddo
1105             do i=nnt,nct
1106               if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1107                 do j=1,3
1108                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1109                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1110                 enddo
1111               endif
1112             enddo
1113             return
1114           else
1115             call read_angles(inp,*36)
1116           endif
1117           goto 37
1118    36     write (iout,'(a)') 'Error reading angle file.'
1119 #ifdef MPI
1120           call mpi_finalize( MPI_COMM_WORLD,IERR )
1121 #endif
1122           stop 'Error reading angle file.'
1123    37     continue
1124         else if (extconf) then
1125          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1126      &    write (iout,'(a)') 'Extended chain initial geometry.'
1127          do i=3,nres
1128           theta(i)=90d0*deg2rad
1129          enddo
1130          do i=4,nres
1131           phi(i)=180d0*deg2rad
1132          enddo
1133          do i=2,nres-1
1134           alph(i)=110d0*deg2rad
1135          enddo
1136          do i=2,nres-1
1137           omeg(i)=-120d0*deg2rad
1138           if (itype(i).le.0) omeg(i)=-omeg(i)
1139          enddo
1140         call chainbuild_extconf
1141         else
1142           if(me.eq.king.or..not.out1file)
1143      &     write (iout,'(a)') 'Random-generated initial geometry.'
1144
1145
1146 #ifdef MPI
1147           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1148      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1149 #endif
1150             do itrial=1,100
1151               itmp=1
1152               call gen_rand_conf(itmp,*30)
1153               goto 40
1154    30         write (iout,*) 'Failed to generate random conformation',
1155      &          ', itrial=',itrial
1156               write (*,*) 'Processor:',me,
1157      &          ' Failed to generate random conformation',
1158      &          ' itrial=',itrial
1159               call intout
1160
1161 #ifdef AIX
1162               call flush_(iout)
1163 #else
1164               call flush(iout)
1165 #endif
1166             enddo
1167             write (iout,'(a,i3,a)') 'Processor:',me,
1168      &        ' error in generating random conformation.'
1169             write (*,'(a,i3,a)') 'Processor:',me,
1170      &        ' error in generating random conformation.'
1171             call flush(iout)
1172 #ifdef MPI
1173             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1174    40       continue
1175           endif
1176 #else
1177           write (*,'(a)')
1178      &      ' error in generating random conformation.'
1179           stop
1180    40     continue
1181 #endif
1182         endif
1183       elseif (modecalc.eq.4) then
1184         read (inp,'(a)') intinname
1185         open (intin,file=intinname,status='old',err=333)
1186         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1187      &  write (iout,'(a)') 'intinname',intinname
1188         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1189         goto 334
1190   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1191 #ifdef MPI
1192         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1193 #endif
1194         stop 'Error opening angle file.'
1195   334   continue
1196
1197       endif
1198 C Generate distance constraints, if the PDB structure is to be regularized.
1199       if (nthread.gt.0) then
1200         call read_threadbase
1201       endif
1202       call setup_var
1203       if (me.eq.king .or. .not. out1file)
1204      & call intout
1205       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1206         write (iout,'(/a,i3,a)')
1207      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1208         write (iout,'(20i4)') (iss(i),i=1,ns)
1209        if (dyn_ss) then
1210           write(iout,*)"Running with dynamic disulfide-bond formation"
1211        else
1212         write (iout,'(/a/)') 'Pre-formed links are:'
1213         do i=1,nss
1214           i1=ihpb(i)-nres
1215           i2=jhpb(i)-nres
1216           it1=itype(i1)
1217           it2=itype(i2)
1218           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1219      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1220      &    ebr,forcon(i)
1221         enddo
1222         write (iout,'(a)')
1223        endif
1224       endif
1225       if (ns.gt.0.and.dyn_ss) then
1226           do i=nss+1,nhpb
1227             ihpb(i-nss)=ihpb(i)
1228             jhpb(i-nss)=jhpb(i)
1229             forcon(i-nss)=forcon(i)
1230             dhpb(i-nss)=dhpb(i)
1231           enddo
1232           nhpb=nhpb-nss
1233           nss=0
1234           call hpb_partition
1235           do i=1,ns
1236             dyn_ss_mask(iss(i))=.true.
1237           enddo
1238       endif
1239       if (i2ndstr.gt.0) call secstrp2dihc
1240 c      call geom_to_var(nvar,x)
1241 c      call etotal(energia(0))
1242 c      call enerprint(energia(0))
1243 c      call briefout(0,etot)
1244 c      stop
1245 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1246 cd    write (iout,'(a)') 'Variable list:'
1247 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1248 #ifdef MPI
1249       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1250      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1251      &  'Processor',myrank,': end reading molecular data.'
1252 #endif
1253       print *,"A TU?"
1254       return
1255       end
1256 c--------------------------------------------------------------------------
1257       logical function seq_comp(itypea,itypeb,length)
1258       implicit none
1259       integer length,itypea(length),itypeb(length)
1260       integer i
1261       do i=1,length
1262         if (itypea(i).ne.itypeb(i)) then
1263           seq_comp=.false.
1264           return
1265         endif
1266       enddo
1267       seq_comp=.true.
1268       return
1269       end
1270 c-----------------------------------------------------------------------------
1271       subroutine read_bridge
1272 C Read information about disulfide bridges.
1273       implicit real*8 (a-h,o-z)
1274       include 'DIMENSIONS'
1275 #ifdef MPI
1276       include 'mpif.h'
1277 #endif
1278       include 'COMMON.IOUNITS'
1279       include 'COMMON.GEO'
1280       include 'COMMON.VAR'
1281       include 'COMMON.INTERACT'
1282       include 'COMMON.LOCAL'
1283       include 'COMMON.NAMES'
1284       include 'COMMON.CHAIN'
1285       include 'COMMON.FFIELD'
1286       include 'COMMON.SBRIDGE'
1287       include 'COMMON.HEADER'
1288       include 'COMMON.CONTROL'
1289       include 'COMMON.DBASE'
1290       include 'COMMON.THREAD'
1291       include 'COMMON.TIME1'
1292       include 'COMMON.SETUP'
1293 C Read bridging residues.
1294       read (inp,*) ns,(iss(i),i=1,ns)
1295       print *,'ns=',ns
1296       if(me.eq.king.or..not.out1file)
1297      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1298 C Check whether the specified bridging residues are cystines.
1299       do i=1,ns
1300         if (itype(iss(i)).ne.1) then
1301           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1302      &   'Do you REALLY think that the residue ',
1303      &    restyp(itype(iss(i))),i,
1304      &   ' can form a disulfide bridge?!!!'
1305           write (*,'(2a,i3,a)')
1306      &   'Do you REALLY think that the residue ',
1307      &    restyp(itype(iss(i))),i,
1308      &   ' can form a disulfide bridge?!!!'
1309 #ifdef MPI
1310          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1311          stop
1312 #endif
1313         endif
1314       enddo
1315 C Read preformed bridges.
1316       if (ns.gt.0) then
1317       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1318       if(fg_rank.eq.0)
1319      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1320       if (nss.gt.0) then
1321         nhpb=nss
1322 C Check if the residues involved in bridges are in the specified list of
1323 C bridging residues.
1324         do i=1,nss
1325           do j=1,i-1
1326             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1327      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1328               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1329      &      ' contains residues present in other pairs.'
1330               write (*,'(a,i3,a)') 'Disulfide pair',i,
1331      &      ' contains residues present in other pairs.'
1332 #ifdef MPI
1333               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1334               stop
1335 #endif
1336             endif
1337           enddo
1338           do j=1,ns
1339             if (ihpb(i).eq.iss(j)) goto 10
1340           enddo
1341           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1342    10     continue
1343           do j=1,ns
1344             if (jhpb(i).eq.iss(j)) goto 20
1345           enddo
1346           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1347    20     continue
1348           dhpb(i)=dbr
1349           forcon(i)=fbr
1350         enddo
1351         do i=1,nss
1352           ihpb(i)=ihpb(i)+nres
1353           jhpb(i)=jhpb(i)+nres
1354         enddo
1355       endif
1356       endif
1357       return
1358       end
1359 c----------------------------------------------------------------------------
1360       subroutine read_x(kanal,*)
1361       implicit real*8 (a-h,o-z)
1362       include 'DIMENSIONS'
1363       include 'COMMON.GEO'
1364       include 'COMMON.VAR'
1365       include 'COMMON.CHAIN'
1366       include 'COMMON.IOUNITS'
1367       include 'COMMON.CONTROL'
1368       include 'COMMON.LOCAL'
1369       include 'COMMON.INTERACT'
1370 c Read coordinates from input
1371 c
1372       read(kanal,'(8f10.5)',end=10,err=10)
1373      &  ((c(l,k),l=1,3),k=1,nres),
1374      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1375       do j=1,3
1376         c(j,nres+1)=c(j,1)
1377         c(j,2*nres)=c(j,nres)
1378       enddo
1379       call int_from_cart1(.false.)
1380       do i=1,nres-1
1381         do j=1,3
1382           dc(j,i)=c(j,i+1)-c(j,i)
1383           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1384         enddo
1385       enddo
1386       do i=nnt,nct
1387         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1388           do j=1,3
1389             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1390             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1391           enddo
1392         endif
1393       enddo
1394
1395       return
1396    10 return1
1397       end
1398 c----------------------------------------------------------------------------
1399       subroutine read_threadbase
1400       implicit real*8 (a-h,o-z)
1401       include 'DIMENSIONS'
1402       include 'COMMON.IOUNITS'
1403       include 'COMMON.GEO'
1404       include 'COMMON.VAR'
1405       include 'COMMON.INTERACT'
1406       include 'COMMON.LOCAL'
1407       include 'COMMON.NAMES'
1408       include 'COMMON.CHAIN'
1409       include 'COMMON.FFIELD'
1410       include 'COMMON.SBRIDGE'
1411       include 'COMMON.HEADER'
1412       include 'COMMON.CONTROL'
1413       include 'COMMON.DBASE'
1414       include 'COMMON.THREAD'
1415       include 'COMMON.TIME1'
1416 C Read pattern database for threading.
1417       read (icbase,*) nseq
1418       do i=1,nseq
1419         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1420      &   nres_base(2,i),nres_base(3,i)
1421         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1422      &   nres_base(1,i))
1423 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1424 c    &   nres_base(2,i),nres_base(3,i)
1425 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1426 c    &   nres_base(1,i))
1427       enddo
1428       close (icbase)
1429       if (weidis.eq.0.0D0) weidis=0.1D0
1430       do i=nnt,nct
1431         do j=i+2,nct
1432           nhpb=nhpb+1
1433           ihpb(nhpb)=i
1434           jhpb(nhpb)=j
1435           forcon(nhpb)=weidis
1436         enddo
1437       enddo
1438       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1439       write (iout,'(a,i5)') 'nexcl: ',nexcl
1440       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1441       return
1442       end
1443 c------------------------------------------------------------------------------
1444       subroutine setup_var
1445       implicit real*8 (a-h,o-z)
1446       include 'DIMENSIONS'
1447       include 'COMMON.IOUNITS'
1448       include 'COMMON.GEO'
1449       include 'COMMON.VAR'
1450       include 'COMMON.INTERACT'
1451       include 'COMMON.LOCAL'
1452       include 'COMMON.NAMES'
1453       include 'COMMON.CHAIN'
1454       include 'COMMON.FFIELD'
1455       include 'COMMON.SBRIDGE'
1456       include 'COMMON.HEADER'
1457       include 'COMMON.CONTROL'
1458       include 'COMMON.DBASE'
1459       include 'COMMON.THREAD'
1460       include 'COMMON.TIME1'
1461 C Set up variable list.
1462       ntheta=nres-2
1463       nphi=nres-3
1464       nvar=ntheta+nphi
1465       nside=0
1466       do i=2,nres-1
1467         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1468           nside=nside+1
1469           ialph(i,1)=nvar+nside
1470           ialph(nside,2)=i
1471         endif
1472       enddo
1473       if (indphi.gt.0) then
1474         nvar=nphi
1475       else if (indback.gt.0) then
1476         nvar=nphi+ntheta
1477       else
1478         nvar=nvar+2*nside
1479       endif
1480 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1481       return
1482       end
1483 c----------------------------------------------------------------------------
1484       subroutine gen_dist_constr
1485 C Generate CA distance constraints.
1486       implicit real*8 (a-h,o-z)
1487       include 'DIMENSIONS'
1488       include 'COMMON.IOUNITS'
1489       include 'COMMON.GEO'
1490       include 'COMMON.VAR'
1491       include 'COMMON.INTERACT'
1492       include 'COMMON.LOCAL'
1493       include 'COMMON.NAMES'
1494       include 'COMMON.CHAIN'
1495       include 'COMMON.FFIELD'
1496       include 'COMMON.SBRIDGE'
1497       include 'COMMON.HEADER'
1498       include 'COMMON.CONTROL'
1499       include 'COMMON.DBASE'
1500       include 'COMMON.THREAD'
1501       include 'COMMON.TIME1'
1502       dimension itype_pdb(maxres)
1503       common /pizda/ itype_pdb
1504       character*2 iden
1505 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1506 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1507 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1508 cd     & ' nsup',nsup
1509       do i=nstart_sup,nstart_sup+nsup-1
1510 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1511 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1512         do j=i+2,nstart_sup+nsup-1
1513           nhpb=nhpb+1
1514           ihpb(nhpb)=i+nstart_seq-nstart_sup
1515           jhpb(nhpb)=j+nstart_seq-nstart_sup
1516           forcon(nhpb)=weidis
1517           dhpb(nhpb)=dist(i,j)
1518         enddo
1519       enddo
1520 cd      write (iout,'(a)') 'Distance constraints:'
1521 cd      do i=nss+1,nhpb
1522 cd        ii=ihpb(i)
1523 cd        jj=jhpb(i)
1524 cd        iden='CA'
1525 cd        if (ii.gt.nres) then
1526 cd          iden='SC'
1527 cd          ii=ii-nres
1528 cd          jj=jj-nres
1529 cd        endif
1530 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1531 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1532 cd     &  dhpb(i),forcon(i)
1533 cd      enddo
1534       return
1535       end
1536 c----------------------------------------------------------------------------
1537       subroutine map_read
1538       implicit real*8 (a-h,o-z)
1539       include 'DIMENSIONS'
1540       include 'COMMON.MAP'
1541       include 'COMMON.IOUNITS'
1542       character*3 angid(4) /'THE','PHI','ALP','OME'/
1543       character*80 mapcard,ucase
1544       do imap=1,nmap
1545         read (inp,'(a)') mapcard
1546         mapcard=ucase(mapcard)
1547         if (index(mapcard,'PHI').gt.0) then
1548           kang(imap)=1
1549         else if (index(mapcard,'THE').gt.0) then
1550           kang(imap)=2
1551         else if (index(mapcard,'ALP').gt.0) then
1552           kang(imap)=3
1553         else if (index(mapcard,'OME').gt.0) then
1554           kang(imap)=4
1555         else
1556           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1557           stop 'Error - illegal variable spec in MAP card.'
1558         endif
1559         call readi (mapcard,'RES1',res1(imap),0)
1560         call readi (mapcard,'RES2',res2(imap),0)
1561         if (res1(imap).eq.0) then
1562           res1(imap)=res2(imap)
1563         else if (res2(imap).eq.0) then
1564           res2(imap)=res1(imap)
1565         endif
1566         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1567           write (iout,'(a)')
1568      &    'Error - illegal definition of variable group in MAP.'
1569           stop 'Error - illegal definition of variable group in MAP.'
1570         endif
1571         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1572         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1573         call readi(mapcard,'NSTEP',nstep(imap),0)
1574         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1575           write (iout,'(a)')
1576      &     'Illegal boundary and/or step size specification in MAP.'
1577           stop 'Illegal boundary and/or step size specification in MAP.'
1578         endif
1579       enddo ! imap
1580       return
1581       end
1582 c----------------------------------------------------------------------------
1583       subroutine csaread
1584       implicit real*8 (a-h,o-z)
1585       include 'DIMENSIONS'
1586       include 'COMMON.IOUNITS'
1587       include 'COMMON.GEO'
1588       include 'COMMON.CSA'
1589       include 'COMMON.BANK'
1590       include 'COMMON.CONTROL'
1591       character*80 ucase
1592       character*620 mcmcard
1593       call card_concat(mcmcard)
1594
1595       call readi(mcmcard,'NCONF',nconf,50)
1596       call readi(mcmcard,'NADD',nadd,0)
1597       call readi(mcmcard,'JSTART',jstart,1)
1598       call readi(mcmcard,'JEND',jend,1)
1599       call readi(mcmcard,'NSTMAX',nstmax,500000)
1600       call readi(mcmcard,'N0',n0,1)
1601       call readi(mcmcard,'N1',n1,6)
1602       call readi(mcmcard,'N2',n2,4)
1603       call readi(mcmcard,'N3',n3,0)
1604       call readi(mcmcard,'N4',n4,0)
1605       call readi(mcmcard,'N5',n5,0)
1606       call readi(mcmcard,'N6',n6,10)
1607       call readi(mcmcard,'N7',n7,0)
1608       call readi(mcmcard,'N8',n8,0)
1609       call readi(mcmcard,'N9',n9,0)
1610       call readi(mcmcard,'N14',n14,0)
1611       call readi(mcmcard,'N15',n15,0)
1612       call readi(mcmcard,'N16',n16,0)
1613       call readi(mcmcard,'N17',n17,0)
1614       call readi(mcmcard,'N18',n18,0)
1615
1616       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1617
1618       call readi(mcmcard,'NDIFF',ndiff,2)
1619       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1620       call readi(mcmcard,'IS1',is1,1)
1621       call readi(mcmcard,'IS2',is2,8)
1622       call readi(mcmcard,'NRAN0',nran0,4)
1623       call readi(mcmcard,'NRAN1',nran1,2)
1624       call readi(mcmcard,'IRR',irr,1)
1625       call readi(mcmcard,'NSEED',nseed,20)
1626       call readi(mcmcard,'NTOTAL',ntotal,10000)
1627       call reada(mcmcard,'CUT1',cut1,2.0d0)
1628       call reada(mcmcard,'CUT2',cut2,5.0d0)
1629       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1630       call readi(mcmcard,'ICMAX',icmax,3)
1631       call readi(mcmcard,'IRESTART',irestart,0)
1632 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1633       ntbankm=0
1634 c!bankt
1635       call reada(mcmcard,'DELE',dele,20.0d0)
1636       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1637       call readi(mcmcard,'IREF',iref,0)
1638       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1639       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1640       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1641       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1642       write (iout,*) "NCONF_IN",nconf_in
1643       return
1644       end
1645 c----------------------------------------------------------------------------
1646 cfmc      subroutine mcmfread
1647 cfmc      implicit real*8 (a-h,o-z)
1648 cfmc      include 'DIMENSIONS'
1649 cfmc      include 'COMMON.MCMF'
1650 cfmc      include 'COMMON.IOUNITS'
1651 cfmc      include 'COMMON.GEO'
1652 cfmc      character*80 ucase
1653 cfmc      character*620 mcmcard
1654 cfmc      call card_concat(mcmcard)
1655 cfmc
1656 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1657 cfmc      write(iout,*)'MAXRANT=',maxrant
1658 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1659 cfmc      write(iout,*)'MAXFAM=',maxfam
1660 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1661 cfmc      write(iout,*)'NNET1=',nnet1
1662 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1663 cfmc      write(iout,*)'NNET2=',nnet2
1664 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1665 cfmc      write(iout,*)'NNET3=',nnet3
1666 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1667 cfmc      write(iout,*)'ILASTT=',ilastt
1668 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1669 cfmc      write(iout,*)'MAXSTR=',maxstr
1670 cfmc      maxstr_f=maxstr/maxfam
1671 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1672 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1673 cfmc      write(iout,*)'NMCMF=',nmcmf
1674 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1675 cfmc      write(iout,*)'IFOCUS=',ifocus
1676 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1677 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1678 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1679 cfmc      write(iout,*)'INTPRT=',intprt
1680 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1681 cfmc      write(iout,*)'IPRT=',iprt
1682 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1683 cfmc      write(iout,*)'IMAXTR=',imaxtr
1684 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1685 cfmc      write(iout,*)'MAXEVEN=',maxeven
1686 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1687 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1688 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1689 cfmc      write(iout,*)'INIMIN=',inimin
1690 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1691 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1692 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1693 cfmc      write(iout,*)'NTHREAD=',nthread
1694 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1695 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1696 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1697 cfmc      write(iout,*)'MAXPERT=',maxpert
1698 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1699 cfmc      write(iout,*)'IRMSD=',irmsd
1700 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1701 cfmc      write(iout,*)'DENEMIN=',denemin
1702 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1703 cfmc      write(iout,*)'RCUT1S=',rcut1s
1704 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1705 cfmc      write(iout,*)'RCUT1E=',rcut1e
1706 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1707 cfmc      write(iout,*)'RCUT2S=',rcut2s
1708 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1709 cfmc      write(iout,*)'RCUT2E=',rcut2e
1710 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1711 cfmc      write(iout,*)'DPERT1=',d_pert1
1712 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1713 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1714 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1715 cfmc      write(iout,*)'DPERT2=',d_pert2
1716 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1717 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1718 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1719 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1720 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1721 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1722 cfmc      d_pert1=deg2rad*d_pert1
1723 cfmc      d_pert1a=deg2rad*d_pert1a
1724 cfmc      d_pert2=deg2rad*d_pert2
1725 cfmc      d_pert2a=deg2rad*d_pert2a
1726 cfmc      d_pert2b=deg2rad*d_pert2b
1727 cfmc      d_pert2c=deg2rad*d_pert2c
1728 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1729 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1730 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1731 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1732 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1733 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1734 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1735 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1736 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1737 cfmc      write(iout,*)'RCUTINI=',rcutini
1738 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1739 cfmc      write(iout,*)'GRAT=',grat
1740 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1741 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1742 cfmc
1743 cfmc      return
1744 cfmc      end
1745 c----------------------------------------------------------------------------
1746       subroutine mcmread
1747       implicit real*8 (a-h,o-z)
1748       include 'DIMENSIONS'
1749       include 'COMMON.MCM'
1750       include 'COMMON.MCE'
1751       include 'COMMON.IOUNITS'
1752       character*80 ucase
1753       character*320 mcmcard
1754       call card_concat(mcmcard)
1755       call readi(mcmcard,'MAXACC',maxacc,100)
1756       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1757       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1758       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1759       call readi(mcmcard,'MAXREPM',maxrepm,200)
1760       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1761       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1762       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1763       call reada(mcmcard,'E_UP',e_up,5.0D0)
1764       call reada(mcmcard,'DELTE',delte,0.1D0)
1765       call readi(mcmcard,'NSWEEP',nsweep,5)
1766       call readi(mcmcard,'NSTEPH',nsteph,0)
1767       call readi(mcmcard,'NSTEPC',nstepc,0)
1768       call reada(mcmcard,'TMIN',tmin,298.0D0)
1769       call reada(mcmcard,'TMAX',tmax,298.0D0)
1770       call readi(mcmcard,'NWINDOW',nwindow,0)
1771       call readi(mcmcard,'PRINT_MC',print_mc,0)
1772       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1773       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1774       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1775       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1776       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1777       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1778       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1779       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1780       if (nwindow.gt.0) then
1781         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1782         do i=1,nwindow
1783           winlen(i)=winend(i)-winstart(i)+1
1784         enddo
1785       endif
1786       if (tmax.lt.tmin) tmax=tmin
1787       if (tmax.eq.tmin) then
1788         nstepc=0
1789         nsteph=0
1790       endif
1791       if (nstepc.gt.0 .and. nsteph.gt.0) then
1792         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1793         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1794       endif
1795 C Probabilities of different move types
1796       sumpro_type(0)=0.0D0
1797       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1798       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1799       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1800       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1801       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1802       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1803       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1804       do i=1,MaxMoveType
1805         print *,'i',i,' sumprotype',sumpro_type(i)
1806         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1807         print *,'i',i,' sumprotype',sumpro_type(i)
1808       enddo
1809       return
1810       end
1811 c----------------------------------------------------------------------------
1812       subroutine read_minim
1813       implicit real*8 (a-h,o-z)
1814       include 'DIMENSIONS'
1815       include 'COMMON.MINIM'
1816       include 'COMMON.IOUNITS'
1817       character*80 ucase
1818       character*320 minimcard
1819       call card_concat(minimcard)
1820       call readi(minimcard,'MAXMIN',maxmin,2000)
1821       call readi(minimcard,'MAXFUN',maxfun,5000)
1822       call readi(minimcard,'MINMIN',minmin,maxmin)
1823       call readi(minimcard,'MINFUN',minfun,maxmin)
1824       call reada(minimcard,'TOLF',tolf,1.0D-2)
1825       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1826       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1827       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1828       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1829       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1830      &         'Options in energy minimization:'
1831       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1832      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1833      & 'MinMin:',MinMin,' MinFun:',MinFun,
1834      & ' TolF:',TolF,' RTolF:',RTolF
1835       return
1836       end
1837 c----------------------------------------------------------------------------
1838       subroutine read_angles(kanal,*)
1839       implicit real*8 (a-h,o-z)
1840       include 'DIMENSIONS'
1841       include 'COMMON.GEO'
1842       include 'COMMON.VAR'
1843       include 'COMMON.CHAIN'
1844       include 'COMMON.IOUNITS'
1845       include 'COMMON.CONTROL'
1846 c Read angles from input
1847 c
1848        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1849        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1850        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1851        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1852
1853        do i=1,nres
1854 c 9/7/01 avoid 180 deg valence angle
1855         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1856 c
1857         theta(i)=deg2rad*theta(i)
1858         phi(i)=deg2rad*phi(i)
1859         alph(i)=deg2rad*alph(i)
1860         omeg(i)=deg2rad*omeg(i)
1861        enddo
1862       return
1863    10 return1
1864       end
1865 c----------------------------------------------------------------------------
1866       subroutine reada(rekord,lancuch,wartosc,default)
1867       implicit none
1868       character*(*) rekord,lancuch
1869       double precision wartosc,default
1870       integer ilen,iread
1871       external ilen
1872       iread=index(rekord,lancuch)
1873       if (iread.eq.0) then
1874         wartosc=default
1875         return
1876       endif
1877       iread=iread+ilen(lancuch)+1
1878       read (rekord(iread:),*,err=10,end=10) wartosc
1879       return
1880   10  wartosc=default
1881       return
1882       end
1883 c----------------------------------------------------------------------------
1884       subroutine readi(rekord,lancuch,wartosc,default)
1885       implicit none
1886       character*(*) rekord,lancuch
1887       integer wartosc,default
1888       integer ilen,iread
1889       external ilen
1890       iread=index(rekord,lancuch)
1891       if (iread.eq.0) then
1892         wartosc=default
1893         return
1894       endif
1895       iread=iread+ilen(lancuch)+1
1896       read (rekord(iread:),*,err=10,end=10) wartosc
1897       return
1898   10  wartosc=default
1899       return
1900       end
1901 c----------------------------------------------------------------------------
1902       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1903       implicit none
1904       integer dim,i
1905       integer tablica(dim),default
1906       character*(*) rekord,lancuch
1907       character*80 aux
1908       integer ilen,iread
1909       external ilen
1910       do i=1,dim
1911         tablica(i)=default
1912       enddo
1913       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1914       if (iread.eq.0) return
1915       iread=iread+ilen(lancuch)+1
1916       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1917    10 return
1918       end
1919 c----------------------------------------------------------------------------
1920       subroutine multreada(rekord,lancuch,tablica,dim,default)
1921       implicit none
1922       integer dim,i
1923       double precision tablica(dim),default
1924       character*(*) rekord,lancuch
1925       character*80 aux
1926       integer ilen,iread
1927       external ilen
1928       do i=1,dim
1929         tablica(i)=default
1930       enddo
1931       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1932       if (iread.eq.0) return
1933       iread=iread+ilen(lancuch)+1
1934       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1935    10 return
1936       end
1937 c----------------------------------------------------------------------------
1938       subroutine openunits
1939       implicit real*8 (a-h,o-z)
1940       include 'DIMENSIONS'
1941 #ifdef MPI
1942       include 'mpif.h'
1943       character*16 form,nodename
1944       integer nodelen
1945 #endif
1946       include 'COMMON.SETUP'
1947       include 'COMMON.IOUNITS'
1948       include 'COMMON.MD'
1949       include 'COMMON.CONTROL'
1950       integer lenpre,lenpot,ilen,lentmp
1951       external ilen
1952       character*3 out1file_text,ucase
1953       character*3 ll
1954       external ucase
1955 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1956       call getenv_loc("PREFIX",prefix)
1957       pref_orig = prefix
1958       call getenv_loc("POT",pot)
1959       call getenv_loc("DIRTMP",tmpdir)
1960       call getenv_loc("CURDIR",curdir)
1961       call getenv_loc("OUT1FILE",out1file_text)
1962 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1963       out1file_text=ucase(out1file_text)
1964       if (out1file_text(1:1).eq."Y") then
1965         out1file=.true.
1966       else
1967         out1file=fg_rank.gt.0
1968       endif
1969       lenpre=ilen(prefix)
1970       lenpot=ilen(pot)
1971       lentmp=ilen(tmpdir)
1972       if (lentmp.gt.0) then
1973           write (*,'(80(1h!))')
1974           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1975           write (*,'(80(1h!))')
1976           write (*,*)"All output files will be on node /tmp directory."
1977 #ifdef MPI
1978         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1979         if (me.eq.king) then
1980           write (*,*) "The master node is ",nodename
1981         else if (fg_rank.eq.0) then
1982           write (*,*) "I am the CG slave node ",nodename
1983         else
1984           write (*,*) "I am the FG slave node ",nodename
1985         endif
1986 #endif
1987         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1988         lenpre = lentmp+lenpre+1
1989       endif
1990       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1991 C Get the names and open the input files
1992 #if defined(WINIFL) || defined(WINPGI)
1993       open(1,file=pref_orig(:ilen(pref_orig))//
1994      &  '.inp',status='old',readonly,shared)
1995        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1996 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1997 C Get parameter filenames and open the parameter files.
1998       call getenv_loc('BONDPAR',bondname)
1999       open (ibond,file=bondname,status='old',readonly,shared)
2000       call getenv_loc('THETPAR',thetname)
2001       open (ithep,file=thetname,status='old',readonly,shared)
2002       call getenv_loc('ROTPAR',rotname)
2003       open (irotam,file=rotname,status='old',readonly,shared)
2004       call getenv_loc('TORPAR',torname)
2005       open (itorp,file=torname,status='old',readonly,shared)
2006       call getenv_loc('TORDPAR',tordname)
2007       open (itordp,file=tordname,status='old',readonly,shared)
2008       call getenv_loc('TORKCC',torkccname)
2009       open (itorkcc,file=torkccname,status='old',readonly,shared)
2010       call getenv_loc('THETKCC',thetkccname)
2011       open (ithetkcc,file=thetkccname,status='old',readonly,shared)
2012       call getenv_loc('FOURIER',fouriername)
2013       open (ifourier,file=fouriername,status='old',readonly,shared)
2014       call getenv_loc('ELEPAR',elename)
2015       open (ielep,file=elename,status='old',readonly,shared)
2016       call getenv_loc('SIDEPAR',sidename)
2017       open (isidep,file=sidename,status='old',readonly,shared)
2018       call getenv_loc('LIPTRANPAR',liptranname)
2019       open (iliptranpar,file=liptranname,status='old',readonly,shared)
2020 #elif (defined CRAY) || (defined AIX)
2021       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2022      &  action='read')
2023 c      print *,"Processor",myrank," opened file 1"
2024       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2025 c      print *,"Processor",myrank," opened file 9"
2026 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2027 C Get parameter filenames and open the parameter files.
2028       call getenv_loc('BONDPAR',bondname)
2029       open (ibond,file=bondname,status='old',action='read')
2030 c      print *,"Processor",myrank," opened file IBOND"
2031       call getenv_loc('THETPAR',thetname)
2032       open (ithep,file=thetname,status='old',action='read')
2033 c      print *,"Processor",myrank," opened file ITHEP"
2034       call getenv_loc('ROTPAR',rotname)
2035       open (irotam,file=rotname,status='old',action='read')
2036 c      print *,"Processor",myrank," opened file IROTAM"
2037       call getenv_loc('TORPAR',torname)
2038       open (itorp,file=torname,status='old',action='read')
2039 c      print *,"Processor",myrank," opened file ITORP"
2040       call getenv_loc('TORDPAR',tordname)
2041       open (itordp,file=tordname,status='old',action='read')
2042       call getenv_loc('TORKCC',torkccname)
2043       open (itorkcc,file=torkccname,status='old',action='read')
2044       call getenv_loc('THETKCC',thetkccname)
2045       open (ithetkcc,file=thetkccname,status='old',action='read')
2046 c      print *,"Processor",myrank," opened file ITORDP"
2047       call getenv_loc('SCCORPAR',sccorname)
2048       open (isccor,file=sccorname,status='old',action='read')
2049 c      print *,"Processor",myrank," opened file ISCCOR"
2050       call getenv_loc('FOURIER',fouriername)
2051       open (ifourier,file=fouriername,status='old',action='read')
2052 c      print *,"Processor",myrank," opened file IFOURIER"
2053       call getenv_loc('ELEPAR',elename)
2054       open (ielep,file=elename,status='old',action='read')
2055 c      print *,"Processor",myrank," opened file IELEP"
2056       call getenv_loc('SIDEPAR',sidename)
2057       open (isidep,file=sidename,status='old',action='read')
2058       call getenv_loc('LIPTRANPAR',liptranname)
2059       open (iliptranpar,file=liptranname,status='old',action='read')
2060 c      print *,"Processor",myrank," opened file ISIDEP"
2061 c      print *,"Processor",myrank," opened parameter files"
2062 #elif (defined G77)
2063       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2064       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2065 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2066 C Get parameter filenames and open the parameter files.
2067       call getenv_loc('BONDPAR',bondname)
2068       open (ibond,file=bondname,status='old')
2069       call getenv_loc('THETPAR',thetname)
2070       open (ithep,file=thetname,status='old')
2071       call getenv_loc('ROTPAR',rotname)
2072       open (irotam,file=rotname,status='old')
2073       call getenv_loc('TORPAR',torname)
2074       open (itorp,file=torname,status='old')
2075       call getenv_loc('TORDPAR',tordname)
2076       open (itordp,file=tordname,status='old')
2077       call getenv_loc('TORKCC',torkccname)
2078       open (itorkcc,file=torkccname,status='old')
2079       call getenv_loc('THETKCC',thetkccname)
2080       open (ithetkcc,file=thetkccname,status='old')
2081       call getenv_loc('SCCORPAR',sccorname)
2082       open (isccor,file=sccorname,status='old')
2083       call getenv_loc('FOURIER',fouriername)
2084       open (ifourier,file=fouriername,status='old')
2085       call getenv_loc('ELEPAR',elename)
2086       open (ielep,file=elename,status='old')
2087       call getenv_loc('SIDEPAR',sidename)
2088       open (isidep,file=sidename,status='old')
2089       call getenv_loc('LIPTRANPAR',liptranname)
2090       open (iliptranpar,file=liptranname,status='old')
2091 #else
2092       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2093      &  readonly)
2094        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2095 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2096 C Get parameter filenames and open the parameter files.
2097       call getenv_loc('BONDPAR',bondname)
2098       open (ibond,file=bondname,status='old',readonly)
2099       call getenv_loc('THETPAR',thetname)
2100       open (ithep,file=thetname,status='old',readonly)
2101       call getenv_loc('ROTPAR',rotname)
2102       open (irotam,file=rotname,status='old',readonly)
2103       call getenv_loc('TORPAR',torname)
2104       open (itorp,file=torname,status='old',readonly)
2105       call getenv_loc('TORDPAR',tordname)
2106       open (itordp,file=tordname,status='old',readonly)
2107       call getenv_loc('TORKCC',torkccname)
2108       open (itorkcc,file=torkccname,status='old',readonly)
2109       call getenv_loc('THETKCC',thetkccname)
2110       open (ithetkcc,file=thetkccname,status='old',readonly)
2111       call getenv_loc('SCCORPAR',sccorname)
2112       open (isccor,file=sccorname,status='old',readonly)
2113 #ifndef CRYST_THETA
2114       call getenv_loc('THETPARPDB',thetname_pdb)
2115       print *,"thetname_pdb ",thetname_pdb
2116       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2117       print *,ithep_pdb," opened"
2118 #endif
2119       call getenv_loc('FOURIER',fouriername)
2120       open (ifourier,file=fouriername,status='old',readonly)
2121       call getenv_loc('ELEPAR',elename)
2122       open (ielep,file=elename,status='old',readonly)
2123       call getenv_loc('SIDEPAR',sidename)
2124       open (isidep,file=sidename,status='old',readonly)
2125       call getenv_loc('LIPTRANPAR',liptranname)
2126       open (iliptranpar,file=liptranname,status='old',action='read')
2127 #ifndef CRYST_SC
2128       call getenv_loc('ROTPARPDB',rotname_pdb)
2129       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2130 #endif
2131 #endif
2132 #ifndef OLDSCP
2133 C
2134 C 8/9/01 In the newest version SCp interaction constants are read from a file
2135 C Use -DOLDSCP to use hard-coded constants instead.
2136 C
2137       call getenv_loc('SCPPAR',scpname)
2138 #if defined(WINIFL) || defined(WINPGI)
2139       open (iscpp,file=scpname,status='old',readonly,shared)
2140 #elif (defined CRAY)  || (defined AIX)
2141       open (iscpp,file=scpname,status='old',action='read')
2142 #elif (defined G77)
2143       open (iscpp,file=scpname,status='old')
2144 #else
2145       open (iscpp,file=scpname,status='old',readonly)
2146 #endif
2147 #endif
2148       call getenv_loc('PATTERN',patname)
2149 #if defined(WINIFL) || defined(WINPGI)
2150       open (icbase,file=patname,status='old',readonly,shared)
2151 #elif (defined CRAY)  || (defined AIX)
2152       open (icbase,file=patname,status='old',action='read')
2153 #elif (defined G77)
2154       open (icbase,file=patname,status='old')
2155 #else
2156       open (icbase,file=patname,status='old',readonly)
2157 #endif
2158 #ifdef MPI
2159 C Open output file only for CG processes
2160 c      print *,"Processor",myrank," fg_rank",fg_rank
2161       if (fg_rank.eq.0) then
2162
2163       if (nodes.eq.1) then
2164         npos=3
2165       else
2166         npos = dlog10(dfloat(nodes-1))+1
2167       endif
2168       if (npos.lt.3) npos=3
2169       write (liczba,'(i1)') npos
2170       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2171      &  //')'
2172       write (liczba,form) me
2173       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2174      &  liczba(:ilen(liczba))
2175       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2176      &  //'.int'
2177       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2178      &  //'.pdb'
2179       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2180      &  liczba(:ilen(liczba))//'.mol2'
2181       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2182      &  liczba(:ilen(liczba))//'.stat'
2183       if (lentmp.gt.0)
2184      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2185      &      //liczba(:ilen(liczba))//'.stat')
2186       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2187      &  //'.rst'
2188       if(usampl) then
2189           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2190      & liczba(:ilen(liczba))//'.const'
2191       endif
2192
2193       endif
2194 #else
2195       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2196       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2197       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2198       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2199       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2200       if (lentmp.gt.0)
2201      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2202      &    //'.stat')
2203       rest2name=prefix(:ilen(prefix))//'.rst'
2204       if(usampl) then
2205          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2206       endif
2207 #endif
2208 #if defined(AIX) || defined(PGI)
2209       if (me.eq.king .or. .not. out1file)
2210      &   open(iout,file=outname,status='unknown')
2211 #ifdef DEBUG
2212       if (fg_rank.gt.0) then
2213         write (liczba,'(i3.3)') myrank/nfgtasks
2214         write (ll,'(bz,i3.3)') fg_rank
2215         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2216      &   status='unknown')
2217       endif
2218 #endif
2219       if(me.eq.king) then
2220        open(igeom,file=intname,status='unknown',position='append')
2221        open(ipdb,file=pdbname,status='unknown')
2222        open(imol2,file=mol2name,status='unknown')
2223        open(istat,file=statname,status='unknown',position='append')
2224       else
2225 c1out       open(iout,file=outname,status='unknown')
2226       endif
2227 #else
2228       if (me.eq.king .or. .not.out1file)
2229      &    open(iout,file=outname,status='unknown')
2230 #ifdef DEBUG
2231       if (fg_rank.gt.0) then
2232         write (liczba,'(i3.3)') myrank/nfgtasks
2233         write (ll,'(bz,i3.3)') fg_rank
2234         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2235      &   status='unknown')
2236       endif
2237 #endif
2238       if(me.eq.king) then
2239        open(igeom,file=intname,status='unknown',access='append')
2240        open(ipdb,file=pdbname,status='unknown')
2241        open(imol2,file=mol2name,status='unknown')
2242        open(istat,file=statname,status='unknown',access='append')
2243       else
2244 c1out       open(iout,file=outname,status='unknown')
2245       endif
2246 #endif
2247       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2248       csa_seed=prefix(:lenpre)//'.CSA.seed'
2249       csa_history=prefix(:lenpre)//'.CSA.history'
2250       csa_bank=prefix(:lenpre)//'.CSA.bank'
2251       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2252       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2253       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2254 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2255       csa_int=prefix(:lenpre)//'.int'
2256       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2257       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2258       csa_in=prefix(:lenpre)//'.CSA.in'
2259 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2260 C Write file names
2261       if (me.eq.king)then
2262       write (iout,'(80(1h-))')
2263       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2264       write (iout,'(80(1h-))')
2265       write (iout,*) "Input file                      : ",
2266      &  pref_orig(:ilen(pref_orig))//'.inp'
2267       write (iout,*) "Output file                     : ",
2268      &  outname(:ilen(outname))
2269       write (iout,*)
2270       write (iout,*) "Sidechain potential file        : ",
2271      &  sidename(:ilen(sidename))
2272 #ifndef OLDSCP
2273       write (iout,*) "SCp potential file              : ",
2274      &  scpname(:ilen(scpname))
2275 #endif
2276       write (iout,*) "Electrostatic potential file    : ",
2277      &  elename(:ilen(elename))
2278       write (iout,*) "Cumulant coefficient file       : ",
2279      &  fouriername(:ilen(fouriername))
2280       write (iout,*) "Torsional parameter file        : ",
2281      &  torname(:ilen(torname))
2282       write (iout,*) "Double torsional parameter file : ",
2283      &  tordname(:ilen(tordname))
2284       write (iout,*) "SCCOR parameter file : ",
2285      &  sccorname(:ilen(sccorname))
2286       write (iout,*) "Bond & inertia constant file    : ",
2287      &  bondname(:ilen(bondname))
2288       write (iout,*) "Bending parameter file          : ",
2289      &  thetname(:ilen(thetname))
2290       write (iout,*) "Rotamer parameter file          : ",
2291      &  rotname(:ilen(rotname))
2292       write (iout,*) "Thetpdb parameter file          : ",
2293      &  thetname_pdb(:ilen(thetname_pdb))
2294       write (iout,*) "Threading database              : ",
2295      &  patname(:ilen(patname))
2296       if (lentmp.ne.0)
2297      &write (iout,*)" DIRTMP                          : ",
2298      &  tmpdir(:lentmp)
2299       write (iout,'(80(1h-))')
2300       endif
2301       return
2302       end
2303 c----------------------------------------------------------------------------
2304       subroutine card_concat(card)
2305       implicit real*8 (a-h,o-z)
2306       include 'DIMENSIONS'
2307       include 'COMMON.IOUNITS'
2308       character*(*) card
2309       character*80 karta,ucase
2310       external ilen
2311       read (inp,'(a)') karta
2312       karta=ucase(karta)
2313       card=' '
2314       do while (karta(80:80).eq.'&')
2315         card=card(:ilen(card)+1)//karta(:79)
2316         read (inp,'(a)') karta
2317         karta=ucase(karta)
2318       enddo
2319       card=card(:ilen(card)+1)//karta
2320       return
2321       end
2322 c----------------------------------------------------------------------------------
2323       subroutine readrst
2324       implicit real*8 (a-h,o-z)
2325       include 'DIMENSIONS'
2326       include 'COMMON.CHAIN'
2327       include 'COMMON.IOUNITS'
2328       include 'COMMON.MD'
2329       open(irest2,file=rest2name,status='unknown')
2330       read(irest2,*) totT,EK,potE,totE,t_bath
2331       totTafm=totT
2332       do i=1,2*nres
2333          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2334       enddo
2335       do i=1,2*nres
2336          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2337       enddo
2338       if(usampl) then
2339              read (irest2,*) iset
2340       endif
2341       close(irest2)
2342       return
2343       end
2344 c---------------------------------------------------------------------------------
2345       subroutine read_fragments
2346       implicit real*8 (a-h,o-z)
2347       include 'DIMENSIONS'
2348 #ifdef MPI
2349       include 'mpif.h'
2350 #endif
2351       include 'COMMON.SETUP'
2352       include 'COMMON.CHAIN'
2353       include 'COMMON.IOUNITS'
2354       include 'COMMON.MD'
2355       include 'COMMON.CONTROL'
2356       read(inp,*) nset,nfrag,npair,nfrag_back
2357       if(me.eq.king.or..not.out1file)
2358      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2359      &  " nfrag_back",nfrag_back
2360       do iset=1,nset
2361          read(inp,*) mset(iset)
2362        do i=1,nfrag
2363          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2364      &     qinfrag(i,iset)
2365          if(me.eq.king.or..not.out1file)
2366      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2367      &     ifrag(2,i,iset), qinfrag(i,iset)
2368        enddo
2369        do i=1,npair
2370         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2371      &    qinpair(i,iset)
2372         if(me.eq.king.or..not.out1file)
2373      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2374      &    ipair(2,i,iset), qinpair(i,iset)
2375        enddo
2376        do i=1,nfrag_back
2377         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2378      &     wfrag_back(3,i,iset),
2379      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2380         if(me.eq.king.or..not.out1file)
2381      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2382      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2383        enddo
2384       enddo
2385       return
2386       end
2387 C---------------------------------------------------------------------------
2388       subroutine read_afminp
2389             implicit real*8 (a-h,o-z)
2390       include 'DIMENSIONS'
2391 #ifdef MPI
2392       include 'mpif.h'
2393 #endif
2394       include 'COMMON.SETUP'
2395       include 'COMMON.CONTROL'
2396       include 'COMMON.CHAIN'
2397       include 'COMMON.IOUNITS'
2398       include 'COMMON.SBRIDGE'
2399       character*320 afmcard
2400       print *, "wchodze"
2401       call card_concat(afmcard)
2402       call readi(afmcard,"BEG",afmbeg,0)
2403       call readi(afmcard,"END",afmend,0)
2404       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2405       call reada(afmcard,"VEL",velAFMconst,0.0d0)
2406       print *,'FORCE=' ,forceAFMconst
2407 CCCC NOW PROPERTIES FOR AFM
2408        distafminit=0.0d0
2409        do i=1,3
2410         distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2411        enddo
2412         distafminit=dsqrt(distafminit)
2413         print *,'initdist',distafminit
2414       return
2415       end
2416
2417 c-------------------------------------------------------------------------------
2418       subroutine read_dist_constr
2419       implicit real*8 (a-h,o-z)
2420       include 'DIMENSIONS'
2421 #ifdef MPI
2422       include 'mpif.h'
2423 #endif
2424       include 'COMMON.SETUP'
2425       include 'COMMON.CONTROL'
2426       include 'COMMON.CHAIN'
2427       include 'COMMON.IOUNITS'
2428       include 'COMMON.SBRIDGE'
2429       integer ifrag_(2,100),ipair_(2,100)
2430       double precision wfrag_(100),wpair_(100)
2431       character*500 controlcard
2432       print *, "WCHODZE"
2433       write (iout,*) "Calling read_dist_constr"
2434 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2435 c      call flush(iout)
2436       call card_concat(controlcard)
2437       call readi(controlcard,"NFRAG",nfrag_,0)
2438       call readi(controlcard,"NPAIR",npair_,0)
2439       call readi(controlcard,"NDIST",ndist_,0)
2440       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2441       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2442       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2443       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2444       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2445 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2446 c      write (iout,*) "IFRAG"
2447 c      do i=1,nfrag_
2448 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2449 c      enddo
2450 c      write (iout,*) "IPAIR"
2451 c      do i=1,npair_
2452 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2453 c      enddo
2454       call flush(iout)
2455       do i=1,nfrag_
2456         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2457         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2458      &    ifrag_(2,i)=nstart_sup+nsup-1
2459 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2460         call flush(iout)
2461         if (wfrag_(i).gt.0.0d0) then
2462         do j=ifrag_(1,i),ifrag_(2,i)-1
2463           do k=j+1,ifrag_(2,i)
2464 c            write (iout,*) "j",j," k",k
2465             ddjk=dist(j,k)
2466             if (constr_dist.eq.1) then
2467             nhpb=nhpb+1
2468             ihpb(nhpb)=j
2469             jhpb(nhpb)=k
2470               dhpb(nhpb)=ddjk
2471             forcon(nhpb)=wfrag_(i)
2472             else if (constr_dist.eq.2) then
2473               if (ddjk.le.dist_cut) then
2474                 nhpb=nhpb+1
2475                 ihpb(nhpb)=j
2476                 jhpb(nhpb)=k
2477                 dhpb(nhpb)=ddjk
2478                 forcon(nhpb)=wfrag_(i)
2479               endif
2480             else
2481               nhpb=nhpb+1
2482               ihpb(nhpb)=j
2483               jhpb(nhpb)=k
2484               dhpb(nhpb)=ddjk
2485               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2486             endif
2487 #ifdef MPI
2488             if (.not.out1file .or. me.eq.king)
2489      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2490      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2491 #else
2492             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2493      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2494 #endif
2495           enddo
2496         enddo
2497         endif
2498       enddo
2499       do i=1,npair_
2500         if (wpair_(i).gt.0.0d0) then
2501         ii = ipair_(1,i)
2502         jj = ipair_(2,i)
2503         if (ii.gt.jj) then
2504           itemp=ii
2505           ii=jj
2506           jj=itemp
2507         endif
2508         do j=ifrag_(1,ii),ifrag_(2,ii)
2509           do k=ifrag_(1,jj),ifrag_(2,jj)
2510             nhpb=nhpb+1
2511             ihpb(nhpb)=j
2512             jhpb(nhpb)=k
2513             forcon(nhpb)=wpair_(i)
2514             dhpb(nhpb)=dist(j,k)
2515 #ifdef MPI
2516             if (.not.out1file .or. me.eq.king)
2517      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2518      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2519 #else
2520             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2521      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2522 #endif
2523           enddo
2524         enddo
2525         endif
2526       enddo
2527       print *,ndist_
2528       do i=1,ndist_
2529         if (constr_dist.eq.11) then
2530         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2531      &     ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2532         fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2533         else
2534 C        print *,"in else"
2535         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2536      &     ibecarb(i),forcon(nhpb+1)
2537         endif
2538         if (forcon(nhpb+1).gt.0.0d0) then
2539           nhpb=nhpb+1
2540           if (ibecarb(i).gt.0) then
2541             ihpb(i)=ihpb(i)+nres
2542             jhpb(i)=jhpb(i)+nres
2543           endif
2544           if (dhpb(nhpb).eq.0.0d0)
2545      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2546         endif
2547 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2548 C        if (forcon(nhpb+1).gt.0.0d0) then
2549 C          nhpb=nhpb+1
2550 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2551 #ifdef MPI
2552           if (.not.out1file .or. me.eq.king)
2553      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2554      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2555 #else
2556           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2557      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2558 #endif
2559
2560       enddo
2561       call flush(iout)
2562       return
2563       end
2564 c-------------------------------------------------------------------------------
2565 #ifdef WINIFL
2566       subroutine flush(iu)
2567       return
2568       end
2569 #endif
2570 #ifdef AIX
2571       subroutine flush(iu)
2572       call flush_(iu)
2573       return
2574       end
2575 #endif
2576 c------------------------------------------------------------------------------
2577       subroutine copy_to_tmp(source)
2578       include "DIMENSIONS"
2579       include "COMMON.IOUNITS"
2580       character*(*) source
2581       character* 256 tmpfile
2582       integer ilen
2583       external ilen
2584       logical ex
2585       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2586       inquire(file=tmpfile,exist=ex)
2587       if (ex) then
2588         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2589      &   " to temporary directory..."
2590         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2591         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2592       endif
2593       return
2594       end
2595 c------------------------------------------------------------------------------
2596       subroutine move_from_tmp(source)
2597       include "DIMENSIONS"
2598       include "COMMON.IOUNITS"
2599       character*(*) source
2600       integer ilen
2601       external ilen
2602       write (*,*) "Moving ",source(:ilen(source)),
2603      & " from temporary directory to working directory"
2604       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2605       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2606       return
2607       end
2608 c------------------------------------------------------------------------------
2609       subroutine random_init(seed)
2610 C
2611 C Initialize random number generator
2612 C
2613       implicit real*8 (a-h,o-z)
2614       include 'DIMENSIONS'
2615 #ifdef MPI
2616       include 'mpif.h'
2617       logical OKRandom, prng_restart
2618       real*8  r1
2619       integer iseed_array(4)
2620 #endif
2621       include 'COMMON.IOUNITS'
2622       include 'COMMON.TIME1'
2623       include 'COMMON.THREAD'
2624       include 'COMMON.SBRIDGE'
2625       include 'COMMON.CONTROL'
2626       include 'COMMON.MCM'
2627       include 'COMMON.MAP'
2628       include 'COMMON.HEADER'
2629       include 'COMMON.CSA'
2630       include 'COMMON.CHAIN'
2631       include 'COMMON.MUCA'
2632       include 'COMMON.MD'
2633       include 'COMMON.FFIELD'
2634       include 'COMMON.SETUP'
2635       iseed=-dint(dabs(seed))
2636       if (iseed.eq.0) then
2637         write (iout,'(/80(1h*)/20x,a/80(1h*))')
2638      &    'Random seed undefined. The program will stop.'
2639         write (*,'(/80(1h*)/20x,a/80(1h*))')
2640      &    'Random seed undefined. The program will stop.'
2641 #ifdef MPI
2642         call mpi_finalize(mpi_comm_world,ierr)
2643 #endif
2644         stop 'Bad random seed.'
2645       endif
2646 #ifdef MPI
2647       if (fg_rank.eq.0) then
2648       seed=seed*(me+1)+1
2649 #ifdef AMD64
2650       if(me.eq.king)
2651      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2652       OKRandom = prng_restart(me,iseed)
2653 #else
2654       do i=1,4
2655        tmp=65536.0d0**(4-i)
2656        iseed_array(i) = dint(seed/tmp)
2657        seed=seed-iseed_array(i)*tmp
2658       enddo
2659       if(me.eq.king)
2660      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2661      &                 (iseed_array(i),i=1,4)
2662       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2663      &                 (iseed_array(i),i=1,4)
2664       OKRandom = prng_restart(me,iseed_array)
2665 #endif
2666       if (OKRandom) then
2667 c        r1 = prng_next(me)
2668         r1=ran_number(0.0D0,1.0D0)
2669         if(me.eq.king)
2670      &   write (iout,*) 'ran_num',r1
2671         if (r1.lt.0.0d0) OKRandom=.false.
2672       endif
2673       if (.not.OKRandom) then
2674         write (iout,*) 'PRNG IS NOT WORKING!!!'
2675         print *,'PRNG IS NOT WORKING!!!'
2676         if (me.eq.0) then
2677          call flush(iout)
2678          call mpi_abort(mpi_comm_world,error_msg,ierr)
2679          stop
2680         else
2681          write (iout,*) 'too many processors for parallel prng'
2682          write (*,*) 'too many processors for parallel prng'
2683          call flush(iout)
2684          stop
2685         endif
2686       endif
2687       endif
2688 #else
2689       call vrndst(iseed)
2690       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
2691 #endif
2692       return
2693       end