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       character*4 sequence(maxres)
 604       integer rescode
 605       double precision x(maxvar)
 606       character*256 pdbfile
 607       character*400 weightcard
 608       character*80 weightcard_t,ucase
 609       dimension itype_pdb(maxres)
 610       common /pizda/ itype_pdb
 611       logical seq_comp,fail
 612       double precision energia(0:n_ene)
 613       integer ilen
 614       external ilen
 615 C
 616 C Body
 617 C
 618 C Read weights of the subsequent energy terms.
 619        call card_concat(weightcard)
 620        call reada(weightcard,'WLONG',wlong,1.0D0)
 621        call reada(weightcard,'WSC',wsc,wlong)
 622        call reada(weightcard,'WSCP',wscp,wlong)
 623        call reada(weightcard,'WELEC',welec,1.0D0)
 624        call reada(weightcard,'WVDWPP',wvdwpp,welec)
 625        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 626        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 627        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 628        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 629        call reada(weightcard,'WTURN3',wturn3,1.0D0)
 630        call reada(weightcard,'WTURN4',wturn4,1.0D0)
 631        call reada(weightcard,'WTURN6',wturn6,1.0D0)
 632        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 633        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 634        call reada(weightcard,'WBOND',wbond,1.0D0)
 635        call reada(weightcard,'WTOR',wtor,1.0D0)
 636        call reada(weightcard,'WTORD',wtor_d,1.0D0)
 637        call reada(weightcard,'WANG',wang,1.0D0)
 638        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 639        call reada(weightcard,'SCAL14',scal14,0.4D0)
 640        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 641        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 642        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 643        call reada(weightcard,'TEMP0',temp0,300.0d0)
 644        call reada(weightcard,'WLT',wliptran,0.0D0)
 645        if (index(weightcard,'SOFT').gt.0) ipot=6
 646 C 12/1/95 Added weight for the multi-body term WCORR
 647        call reada(weightcard,'WCORRH',wcorr,1.0D0)
 648        if (wcorr4.gt.0.0d0) wcorr=wcorr4
 649        weights(1)=wsc
 650        weights(2)=wscp
 651        weights(3)=welec
 652        weights(4)=wcorr
 653        weights(5)=wcorr5
 654        weights(6)=wcorr6
 655        weights(7)=wel_loc
 656        weights(8)=wturn3
 657        weights(9)=wturn4
 658        weights(10)=wturn6
 659        weights(11)=wang
 660        weights(12)=wscloc
 661        weights(13)=wtor
 662        weights(14)=wtor_d
 663        weights(15)=wstrain
 664        weights(16)=wvdwpp
 665        weights(17)=wbond
 666        weights(18)=scal14
 667        weights(21)=wsccor
 668       if(me.eq.king.or..not.out1file)
 669      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 670      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 671      &  wturn4,wturn6
 672    10 format (/'Energy-term weights (unscaled):'//
 673      & 'WSCC=   ',f10.6,' (SC-SC)'/
 674      & 'WSCP=   ',f10.6,' (SC-p)'/
 675      & 'WELEC=  ',f10.6,' (p-p electr)'/
 676      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 677      & 'WBOND=  ',f10.6,' (stretching)'/
 678      & 'WANG=   ',f10.6,' (bending)'/
 679      & 'WSCLOC= ',f10.6,' (SC local)'/
 680      & 'WTOR=   ',f10.6,' (torsional)'/
 681      & 'WTORD=  ',f10.6,' (double torsional)'/
 682      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 683      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 684      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 685      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 686      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 687      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
 688      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 689      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 690      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 691       if(me.eq.king.or..not.out1file)then
 692        if (wcorr4.gt.0.0d0) then
 693         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 694      &   'between contact pairs of peptide groups'
 695         write (iout,'(2(a,f5.3/))')
 696      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 697      &  'Range of quenching the correlation terms:',2*delt_corr
 698        else if (wcorr.gt.0.0d0) then
 699         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 700      &   'between contact pairs of peptide groups'
 701        endif
 702        write (iout,'(a,f8.3)')
 703      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 704        write (iout,'(a,f8.3)')
 705      &  'General scaling factor of SC-p interactions:',scalscp
 706       endif
 707       r0_corr=cutoff_corr-delt_corr
 708       do i=1,ntyp
 709         aad(i,1)=scalscp*aad(i,1)
 710         aad(i,2)=scalscp*aad(i,2)
 711         bad(i,1)=scalscp*bad(i,1)
 712         bad(i,2)=scalscp*bad(i,2)
 713       enddo
 714       call rescale_weights(t_bath)
 715       if(me.eq.king.or..not.out1file)
 716      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 717      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 718      &  wturn4,wturn6
 719    22 format (/'Energy-term weights (scaled):'//
 720      & 'WSCC=   ',f10.6,' (SC-SC)'/
 721      & 'WSCP=   ',f10.6,' (SC-p)'/
 722      & 'WELEC=  ',f10.6,' (p-p electr)'/
 723      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 724      & 'WBOND=  ',f10.6,' (stretching)'/
 725      & 'WANG=   ',f10.6,' (bending)'/
 726      & 'WSCLOC= ',f10.6,' (SC local)'/
 727      & 'WTOR=   ',f10.6,' (torsional)'/
 728      & 'WTORD=  ',f10.6,' (double torsional)'/
 729      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 730      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 731      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 732      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 733      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 734      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
 735      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 736      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 737      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 738       if(me.eq.king.or..not.out1file)
 739      & write (iout,*) "Reference temperature for weights calculation:",
 740      &  temp0
 741       call reada(weightcard,"D0CM",d0cm,3.78d0)
 742       call reada(weightcard,"AKCM",akcm,15.1d0)
 743       call reada(weightcard,"AKTH",akth,11.0d0)
 744       call reada(weightcard,"AKCT",akct,12.0d0)
 745       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 746       call reada(weightcard,"V2SS",v2ss,7.61d0)
 747       call reada(weightcard,"V3SS",v3ss,13.7d0)
 748       call reada(weightcard,"EBR",ebr,-5.50D0)
 749       call reada(weightcard,"ATRISS",atriss,0.301D0)
 750       call reada(weightcard,"BTRISS",btriss,0.021D0)
 751       call reada(weightcard,"CTRISS",ctriss,1.001D0)
 752       call reada(weightcard,"DTRISS",dtriss,1.001D0)
 753       write (iout,*) "ATRISS=", atriss
 754       write (iout,*) "BTRISS=", btriss
 755       write (iout,*) "CTRISS=", ctriss
 756       write (iout,*) "DTRISS=", dtriss
 757       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 758       do i=1,maxres
 759         dyn_ss_mask(i)=.false.
 760       enddo
 761       do i=1,maxres-1
 762         do j=i+1,maxres
 763           dyn_ssbond_ij(i,j)=1.0d300
 764         enddo
 765       enddo
 766       call reada(weightcard,"HT",Ht,0.0D0)
 767       if (dyn_ss) then
 768         ss_depth=ebr/wsc-0.25*eps(1,1)
 769         Ht=Ht/wsc-0.25*eps(1,1)
 770         akcm=akcm*wstrain/wsc
 771         akth=akth*wstrain/wsc
 772         akct=akct*wstrain/wsc
 773         v1ss=v1ss*wstrain/wsc
 774         v2ss=v2ss*wstrain/wsc
 775         v3ss=v3ss*wstrain/wsc
 776       else
 777         if (wstrain.ne.0.0) then
 778          ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 779         else
 780           ss_depth=0.0
 781         endif
 782       endif
 783
 784       if(me.eq.king.or..not.out1file) then
 785        write (iout,*) "Parameters of the SS-bond potential:"
 786        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 787      & " AKCT",akct
 788        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 789        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 790        write (iout,*)" HT",Ht
 791        print *,'indpdb=',indpdb,' pdbref=',pdbref
 792       endif
 793       if (indpdb.gt.0 .or. pdbref) then
 794         read(inp,'(a)') pdbfile
 795         if(me.eq.king.or..not.out1file)
 796      &   write (iout,'(2a)') 'PDB data will be read from file ',
 797      &   pdbfile(:ilen(pdbfile))
 798         open(ipdbin,file=pdbfile,status='old',err=33)
 799         goto 34
 800   33    write (iout,'(a)') 'Error opening PDB file.'
 801         stop
 802   34    continue
 803 c        write (iout,*) 'Begin reading pdb data'
 804 c        call flush(iout)
 805         call readpdb
 806 c        write (iout,*) 'Finished reading pdb data'
 807 c        call flush(iout)
 808         if(me.eq.king.or..not.out1file)
 809      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 810      &   ' nstart_sup=',nstart_sup
 811         do i=1,nres
 812           itype_pdb(i)=itype(i)
 813         enddo
 814         close (ipdbin)
 815         nnt=nstart_sup
 816         nct=nstart_sup+nsup-1
 817         call contact(.false.,ncont_ref,icont_ref,co)
 818
 819         if (sideadd) then
 820          if(me.eq.king.or..not.out1file)
 821      &    write(iout,*)'Adding sidechains'
 822          maxsi=1000
 823          do i=2,nres-1
 824           iti=itype(i)
 825           if (iti.ne.10 .and. itype(i).ne.ntyp1) then
 826             nsi=0
 827             fail=.true.
 828             do while (fail.and.nsi.le.maxsi)
 829               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 830               nsi=nsi+1
 831             enddo
 832             if(fail) write(iout,*)'Adding sidechain failed for res ',
 833      &              i,' after ',nsi,' trials'
 834           endif
 835          enddo
 836         endif
 837       endif
 838       if (indpdb.eq.0) then
 839 C Read sequence if not taken from the pdb file.
 840         read (inp,*) nres
 841 c        print *,'nres=',nres
 842         if (iscode.gt.0) then
 843           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 844         else
 845           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 846         endif
 847 C Convert sequence to numeric code
 848         do i=1,nres
 849           itype(i)=rescode(i,sequence(i),iscode)
 850         enddo
 851 C Assign initial virtual bond lengths
 852         do i=2,nres
 853           vbld(i)=vbl
 854           vbld_inv(i)=vblinv
 855         enddo
 856         do i=2,nres-1
 857           vbld(i+nres)=dsc(iabs(itype(i)))
 858           vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
 859 c          write (iout,*) "i",i," itype",itype(i),
 860 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 861         enddo
 862       endif
 863 c      print *,nres
 864 c      print '(20i4)',(itype(i),i=1,nres)
 865       do i=1,nres
 866 #ifdef PROCOR
 867         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 868 #else
 869         if (itype(i).eq.ntyp1) then
 870 #endif
 871           itel(i)=0
 872 #ifdef PROCOR
 873         else if (iabs(itype(i+1)).ne.20) then
 874 #else
 875         else if (iabs(itype(i)).ne.20) then
 876 #endif
 877           itel(i)=1
 878         else
 879           itel(i)=2
 880         endif
 881       enddo
 882       if(me.eq.king.or..not.out1file)then
 883        write (iout,*) "ITEL"
 884        do i=1,nres-1
 885          write (iout,*) i,itype(i),itel(i)
 886        enddo
 887        print *,'Call Read_Bridge.'
 888       endif
 889       call read_bridge
 890 C 8/13/98 Set limits to generating the dihedral angles
 891       do i=1,nres
 892         phibound(1,i)=-pi
 893         phibound(2,i)=pi
 894       enddo
 895       read (inp,*) ndih_constr
 896       if (ndih_constr.gt.0) then
 897 C        read (inp,*) ftors
 898         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
 899      &  i=1,ndih_constr)
 900         if(me.eq.king.or..not.out1file)then
 901          write (iout,*)
 902      &   'There are',ndih_constr,' constraints on phi angles.'
 903          do i=1,ndih_constr
 904           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
 905      &    ftors(i)
 906          enddo
 907         endif
 908         do i=1,ndih_constr
 909           phi0(i)=deg2rad*phi0(i)
 910           drange(i)=deg2rad*drange(i)
 911         enddo
 912 C        if(me.eq.king.or..not.out1file)
 913 C     &   write (iout,*) 'FTORS',ftors
 914         do i=1,ndih_constr
 915           ii = idih_constr(i)
 916           phibound(1,ii) = phi0(i)-drange(i)
 917           phibound(2,ii) = phi0(i)+drange(i)
 918         enddo
 919       endif
 920 C first setting the theta boundaries to 0 to pi
 921 C this mean that there is no energy penalty for any angle occuring this can be applied
 922 C for generate random conformation but is not implemented in this way
 923 C      do i=1,nres
 924 C        thetabound(1,i)=0
 925 C        thetabound(2,i)=pi
 926 C      enddo
 927 C begin reading theta constrains this is quartic constrains allowing to
 928 C have smooth second derivative
 929       if (with_theta_constr) then
 930 C with_theta_constr is keyword allowing for occurance of theta constrains
 931       read (inp,*) ntheta_constr
 932 C ntheta_constr is the number of theta constrains
 933       if (ntheta_constr.gt.0) then
 934 C        read (inp,*) ftors
 935         read (inp,*) (itheta_constr(i),theta_constr0(i),
 936      &  theta_drange(i),for_thet_constr(i),
 937      &  i=1,ntheta_constr)
 938 C the above code reads from 1 to ntheta_constr
 939 C itheta_constr(i) residue i for which is theta_constr
 940 C theta_constr0 the global minimum value
 941 C theta_drange is range for which there is no energy penalty
 942 C for_thet_constr is the force constant for quartic energy penalty
 943 C E=k*x**4
 944         if(me.eq.king.or..not.out1file)then
 945          write (iout,*)
 946      &   'There are',ntheta_constr,' constraints on phi angles.'
 947          do i=1,ntheta_constr
 948           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
 949      &    theta_drange(i),
 950      &    for_thet_constr(i)
 951          enddo
 952         endif
 953         do i=1,ntheta_constr
 954           theta_constr0(i)=deg2rad*theta_constr0(i)
 955           theta_drange(i)=deg2rad*theta_drange(i)
 956         enddo
 957 C        if(me.eq.king.or..not.out1file)
 958 C     &   write (iout,*) 'FTORS',ftors
 959 C        do i=1,ntheta_constr
 960 C          ii = itheta_constr(i)
 961 C          thetabound(1,ii) = phi0(i)-drange(i)
 962 C          thetabound(2,ii) = phi0(i)+drange(i)
 963 C        enddo
 964       endif ! ntheta_constr.gt.0
 965       endif! with_theta_constr
 966 C
 967 C      with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
 968 C      write (iout,*) "with_dihed_constr ",with_dihed_constr
 969       nnt=1
 970 #ifdef MPI
 971       if (me.eq.king) then
 972 #endif
 973        write (iout,'(a)') 'Boundaries in phi angle sampling:'
 974        do i=1,nres
 975          write (iout,'(a3,i5,2f10.1)')
 976      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 977        enddo
 978 #ifdef MP
 979       endif
 980 #endif
 981       nct=nres
 982 cd      print *,'NNT=',NNT,' NCT=',NCT
 983       if (itype(1).eq.ntyp1) nnt=2
 984       if (itype(nres).eq.ntyp1) nct=nct-1
 985       if (pdbref) then
 986         if(me.eq.king.or..not.out1file)
 987      &   write (iout,'(a,i3)') 'nsup=',nsup
 988         nstart_seq=nnt
 989         if (nsup.le.(nct-nnt+1)) then
 990           do i=0,nct-nnt+1-nsup
 991             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 992               nstart_seq=nnt+i
 993               goto 111
 994             endif
 995           enddo
 996           write (iout,'(a)')
 997      &            'Error - sequences to be superposed do not match.'
 998           stop
 999         else
1000           do i=0,nsup-(nct-nnt+1)
1001             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1002      &      then
1003               nstart_sup=nstart_sup+i
1004               nsup=nct-nnt+1
1005               goto 111
1006             endif
1007           enddo
1008           write (iout,'(a)')
1009      &            'Error - sequences to be superposed do not match.'
1010         endif
1011   111   continue
1012         if (nsup.eq.0) nsup=nct-nnt
1013         if (nstart_sup.eq.0) nstart_sup=nnt
1014         if (nstart_seq.eq.0) nstart_seq=nnt
1015         if(me.eq.king.or..not.out1file)
1016      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1017      &                 ' nstart_seq=',nstart_seq
1018       endif
1019 c--- Zscore rms -------
1020       if (nz_start.eq.0) nz_start=nnt
1021       if (nz_end.eq.0 .and. nsup.gt.0) then
1022         nz_end=nnt+nsup-1
1023       else if (nz_end.eq.0) then
1024         nz_end=nct
1025       endif
1026       if(me.eq.king.or..not.out1file)then
1027        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1028        write (iout,*) 'IZ_SC=',iz_sc
1029       endif
1030 c----------------------
1031       call init_int_table
1032       if (refstr) then
1033         if (.not.pdbref) then
1034           call read_angles(inp,*38)
1035           goto 39
1036    38     write (iout,'(a)') 'Error reading reference structure.'
1037 #ifdef MPI
1038           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1039           stop 'Error reading reference structure'
1040 #endif
1041    39     call chainbuild
1042           call setup_var
1043 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1044           nstart_sup=nnt
1045           nstart_seq=nnt
1046           nsup=nct-nnt+1
1047           kkk=1
1048           do i=1,2*nres
1049             do j=1,3
1050               cref(j,i,kkk)=c(j,i)
1051             enddo
1052           enddo
1053           call contact(.true.,ncont_ref,icont_ref,co)
1054         endif
1055         endif
1056         print *, "A TU"
1057         write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1058         call flush(iout)
1059         if (constr_dist.gt.0) call read_dist_constr
1060         write (iout,*) "After read_dist_constr nhpb",nhpb
1061         if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1062         call hpb_partition
1063         if(me.eq.king.or..not.out1file)
1064      &   write (iout,*) 'Contact order:',co
1065         if (pdbref) then
1066         if(me.eq.king.or..not.out1file)
1067      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1068         do i=1,ncont_ref
1069           do j=1,2
1070             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1071           enddo
1072           if(me.eq.king.or..not.out1file)
1073      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1074      &     icont_ref(1,i),' ',
1075      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1076         enddo
1077         endif
1078 C      endif
1079       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1080      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1081      &    modecalc.ne.10) then
1082 C If input structure hasn't been supplied from the PDB file read or generate
1083 C initial geometry.
1084         if (iranconf.eq.0 .and. .not. extconf) then
1085           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1086      &     write (iout,'(a)') 'Initial geometry will be read in.'
1087           if (read_cart) then
1088             read(inp,'(8f10.5)',end=36,err=36)
1089      &       ((c(l,k),l=1,3),k=1,nres),
1090      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1091             write (iout,*) "Exit READ_CART"
1092             write (iout,'(8f10.5)')
1093      &       ((c(l,k),l=1,3),k=1,nres),
1094      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1095             call int_from_cart1(.true.)
1096             write (iout,*) "Finish INT_TO_CART"
1097             do i=1,nres-1
1098               do j=1,3
1099                 dc(j,i)=c(j,i+1)-c(j,i)
1100                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1101               enddo
1102             enddo
1103             do i=nnt,nct
1104               if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1105                 do j=1,3
1106                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1107                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1108                 enddo
1109               endif
1110             enddo
1111             return
1112           else
1113             call read_angles(inp,*36)
1114           endif
1115           goto 37
1116    36     write (iout,'(a)') 'Error reading angle file.'
1117 #ifdef MPI
1118           call mpi_finalize( MPI_COMM_WORLD,IERR )
1119 #endif
1120           stop 'Error reading angle file.'
1121    37     continue
1122         else if (extconf) then
1123          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1124      &    write (iout,'(a)') 'Extended chain initial geometry.'
1125          do i=3,nres
1126           theta(i)=90d0*deg2rad
1127          enddo
1128          do i=4,nres
1129           phi(i)=180d0*deg2rad
1130          enddo
1131          do i=2,nres-1
1132           alph(i)=110d0*deg2rad
1133          enddo
1134          do i=2,nres-1
1135           omeg(i)=-120d0*deg2rad
1136           if (itype(i).le.0) omeg(i)=-omeg(i)
1137          enddo
1138         call chainbuild_extconf
1139         else
1140           if(me.eq.king.or..not.out1file)
1141      &     write (iout,'(a)') 'Random-generated initial geometry.'
1142
1143
1144 #ifdef MPI
1145           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1146      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1147 #endif
1148             do itrial=1,100
1149               itmp=1
1150               call gen_rand_conf(itmp,*30)
1151               goto 40
1152    30         write (iout,*) 'Failed to generate random conformation',
1153      &          ', itrial=',itrial
1154               write (*,*) 'Processor:',me,
1155      &          ' Failed to generate random conformation',
1156      &          ' itrial=',itrial
1157               call intout
1158
1159 #ifdef AIX
1160               call flush_(iout)
1161 #else
1162               call flush(iout)
1163 #endif
1164             enddo
1165             write (iout,'(a,i3,a)') 'Processor:',me,
1166      &        ' error in generating random conformation.'
1167             write (*,'(a,i3,a)') 'Processor:',me,
1168      &        ' error in generating random conformation.'
1169             call flush(iout)
1170 #ifdef MPI
1171             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1172    40       continue
1173           endif
1174 #else
1175           write (*,'(a)')
1176      &      ' error in generating random conformation.'
1177           stop
1178    40     continue
1179 #endif
1180         endif
1181       elseif (modecalc.eq.4) then
1182         read (inp,'(a)') intinname
1183         open (intin,file=intinname,status='old',err=333)
1184         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1185      &  write (iout,'(a)') 'intinname',intinname
1186         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1187         goto 334
1188   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1189 #ifdef MPI
1190         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1191 #endif
1192         stop 'Error opening angle file.'
1193   334   continue
1194
1195       endif
1196 C Generate distance constraints, if the PDB structure is to be regularized.
1197       if (nthread.gt.0) then
1198         call read_threadbase
1199       endif
1200       call setup_var
1201       if (me.eq.king .or. .not. out1file)
1202      & call intout
1203       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1204         write (iout,'(/a,i3,a)')
1205      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1206         write (iout,'(20i4)') (iss(i),i=1,ns)
1207        if (dyn_ss) then
1208           write(iout,*)"Running with dynamic disulfide-bond formation"
1209        else
1210         write (iout,'(/a/)') 'Pre-formed links are:'
1211         do i=1,nss
1212           i1=ihpb(i)-nres
1213           i2=jhpb(i)-nres
1214           it1=itype(i1)
1215           it2=itype(i2)
1216           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1217      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1218      &    ebr,forcon(i)
1219         enddo
1220         write (iout,'(a)')
1221        endif
1222       endif
1223       if (ns.gt.0.and.dyn_ss) then
1224           do i=nss+1,nhpb
1225             ihpb(i-nss)=ihpb(i)
1226             jhpb(i-nss)=jhpb(i)
1227             forcon(i-nss)=forcon(i)
1228             dhpb(i-nss)=dhpb(i)
1229           enddo
1230           nhpb=nhpb-nss
1231           nss=0
1232           call hpb_partition
1233           do i=1,ns
1234             dyn_ss_mask(iss(i))=.true.
1235           enddo
1236       endif
1237       if (i2ndstr.gt.0) call secstrp2dihc
1238 c      call geom_to_var(nvar,x)
1239 c      call etotal(energia(0))
1240 c      call enerprint(energia(0))
1241 c      call briefout(0,etot)
1242 c      stop
1243 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1244 cd    write (iout,'(a)') 'Variable list:'
1245 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1246 #ifdef MPI
1247       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1248      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1249      &  'Processor',myrank,': end reading molecular data.'
1250 #endif
1251       print *,"A TU?"
1252       return
1253       end
1254 c--------------------------------------------------------------------------
1255       logical function seq_comp(itypea,itypeb,length)
1256       implicit none
1257       integer length,itypea(length),itypeb(length)
1258       integer i
1259       do i=1,length
1260         if (itypea(i).ne.itypeb(i)) then
1261           seq_comp=.false.
1262           return
1263         endif
1264       enddo
1265       seq_comp=.true.
1266       return
1267       end
1268 c-----------------------------------------------------------------------------
1269       subroutine read_bridge
1270 C Read information about disulfide bridges.
1271       implicit real*8 (a-h,o-z)
1272       include 'DIMENSIONS'
1273 #ifdef MPI
1274       include 'mpif.h'
1275 #endif
1276       include 'COMMON.IOUNITS'
1277       include 'COMMON.GEO'
1278       include 'COMMON.VAR'
1279       include 'COMMON.INTERACT'
1280       include 'COMMON.LOCAL'
1281       include 'COMMON.NAMES'
1282       include 'COMMON.CHAIN'
1283       include 'COMMON.FFIELD'
1284       include 'COMMON.SBRIDGE'
1285       include 'COMMON.HEADER'
1286       include 'COMMON.CONTROL'
1287       include 'COMMON.DBASE'
1288       include 'COMMON.THREAD'
1289       include 'COMMON.TIME1'
1290       include 'COMMON.SETUP'
1291 C Read bridging residues.
1292       read (inp,*) ns,(iss(i),i=1,ns)
1293       print *,'ns=',ns
1294       if(me.eq.king.or..not.out1file)
1295      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1296 C Check whether the specified bridging residues are cystines.
1297       do i=1,ns
1298         if (itype(iss(i)).ne.1) then
1299           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1300      &   'Do you REALLY think that the residue ',
1301      &    restyp(itype(iss(i))),i,
1302      &   ' can form a disulfide bridge?!!!'
1303           write (*,'(2a,i3,a)')
1304      &   'Do you REALLY think that the residue ',
1305      &    restyp(itype(iss(i))),i,
1306      &   ' can form a disulfide bridge?!!!'
1307 #ifdef MPI
1308          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1309          stop
1310 #endif
1311         endif
1312       enddo
1313 C Read preformed bridges.
1314       if (ns.gt.0) then
1315       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1316       if(fg_rank.eq.0)
1317      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1318       if (nss.gt.0) then
1319         nhpb=nss
1320 C Check if the residues involved in bridges are in the specified list of
1321 C bridging residues.
1322         do i=1,nss
1323           do j=1,i-1
1324             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1325      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1326               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1327      &      ' contains residues present in other pairs.'
1328               write (*,'(a,i3,a)') 'Disulfide pair',i,
1329      &      ' contains residues present in other pairs.'
1330 #ifdef MPI
1331               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1332               stop
1333 #endif
1334             endif
1335           enddo
1336           do j=1,ns
1337             if (ihpb(i).eq.iss(j)) goto 10
1338           enddo
1339           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1340    10     continue
1341           do j=1,ns
1342             if (jhpb(i).eq.iss(j)) goto 20
1343           enddo
1344           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1345    20     continue
1346           dhpb(i)=dbr
1347           forcon(i)=fbr
1348         enddo
1349         do i=1,nss
1350           ihpb(i)=ihpb(i)+nres
1351           jhpb(i)=jhpb(i)+nres
1352         enddo
1353       endif
1354       endif
1355       return
1356       end
1357 c----------------------------------------------------------------------------
1358       subroutine read_x(kanal,*)
1359       implicit real*8 (a-h,o-z)
1360       include 'DIMENSIONS'
1361       include 'COMMON.GEO'
1362       include 'COMMON.VAR'
1363       include 'COMMON.CHAIN'
1364       include 'COMMON.IOUNITS'
1365       include 'COMMON.CONTROL'
1366       include 'COMMON.LOCAL'
1367       include 'COMMON.INTERACT'
1368 c Read coordinates from input
1369 c
1370       read(kanal,'(8f10.5)',end=10,err=10)
1371      &  ((c(l,k),l=1,3),k=1,nres),
1372      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1373       do j=1,3
1374         c(j,nres+1)=c(j,1)
1375         c(j,2*nres)=c(j,nres)
1376       enddo
1377       call int_from_cart1(.false.)
1378       do i=1,nres-1
1379         do j=1,3
1380           dc(j,i)=c(j,i+1)-c(j,i)
1381           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1382         enddo
1383       enddo
1384       do i=nnt,nct
1385         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1386           do j=1,3
1387             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1388             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1389           enddo
1390         endif
1391       enddo
1392
1393       return
1394    10 return1
1395       end
1396 c----------------------------------------------------------------------------
1397       subroutine read_threadbase
1398       implicit real*8 (a-h,o-z)
1399       include 'DIMENSIONS'
1400       include 'COMMON.IOUNITS'
1401       include 'COMMON.GEO'
1402       include 'COMMON.VAR'
1403       include 'COMMON.INTERACT'
1404       include 'COMMON.LOCAL'
1405       include 'COMMON.NAMES'
1406       include 'COMMON.CHAIN'
1407       include 'COMMON.FFIELD'
1408       include 'COMMON.SBRIDGE'
1409       include 'COMMON.HEADER'
1410       include 'COMMON.CONTROL'
1411       include 'COMMON.DBASE'
1412       include 'COMMON.THREAD'
1413       include 'COMMON.TIME1'
1414 C Read pattern database for threading.
1415       read (icbase,*) nseq
1416       do i=1,nseq
1417         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1418      &   nres_base(2,i),nres_base(3,i)
1419         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1420      &   nres_base(1,i))
1421 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1422 c    &   nres_base(2,i),nres_base(3,i)
1423 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1424 c    &   nres_base(1,i))
1425       enddo
1426       close (icbase)
1427       if (weidis.eq.0.0D0) weidis=0.1D0
1428       do i=nnt,nct
1429         do j=i+2,nct
1430           nhpb=nhpb+1
1431           ihpb(nhpb)=i
1432           jhpb(nhpb)=j
1433           forcon(nhpb)=weidis
1434         enddo
1435       enddo
1436       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1437       write (iout,'(a,i5)') 'nexcl: ',nexcl
1438       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1439       return
1440       end
1441 c------------------------------------------------------------------------------
1442       subroutine setup_var
1443       implicit real*8 (a-h,o-z)
1444       include 'DIMENSIONS'
1445       include 'COMMON.IOUNITS'
1446       include 'COMMON.GEO'
1447       include 'COMMON.VAR'
1448       include 'COMMON.INTERACT'
1449       include 'COMMON.LOCAL'
1450       include 'COMMON.NAMES'
1451       include 'COMMON.CHAIN'
1452       include 'COMMON.FFIELD'
1453       include 'COMMON.SBRIDGE'
1454       include 'COMMON.HEADER'
1455       include 'COMMON.CONTROL'
1456       include 'COMMON.DBASE'
1457       include 'COMMON.THREAD'
1458       include 'COMMON.TIME1'
1459 C Set up variable list.
1460       ntheta=nres-2
1461       nphi=nres-3
1462       nvar=ntheta+nphi
1463       nside=0
1464       do i=2,nres-1
1465         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1466           nside=nside+1
1467           ialph(i,1)=nvar+nside
1468           ialph(nside,2)=i
1469         endif
1470       enddo
1471       if (indphi.gt.0) then
1472         nvar=nphi
1473       else if (indback.gt.0) then
1474         nvar=nphi+ntheta
1475       else
1476         nvar=nvar+2*nside
1477       endif
1478 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1479       return
1480       end
1481 c----------------------------------------------------------------------------
1482       subroutine gen_dist_constr
1483 C Generate CA distance constraints.
1484       implicit real*8 (a-h,o-z)
1485       include 'DIMENSIONS'
1486       include 'COMMON.IOUNITS'
1487       include 'COMMON.GEO'
1488       include 'COMMON.VAR'
1489       include 'COMMON.INTERACT'
1490       include 'COMMON.LOCAL'
1491       include 'COMMON.NAMES'
1492       include 'COMMON.CHAIN'
1493       include 'COMMON.FFIELD'
1494       include 'COMMON.SBRIDGE'
1495       include 'COMMON.HEADER'
1496       include 'COMMON.CONTROL'
1497       include 'COMMON.DBASE'
1498       include 'COMMON.THREAD'
1499       include 'COMMON.TIME1'
1500       dimension itype_pdb(maxres)
1501       common /pizda/ itype_pdb
1502       character*2 iden
1503 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1504 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1505 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1506 cd     & ' nsup',nsup
1507       do i=nstart_sup,nstart_sup+nsup-1
1508 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1509 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1510         do j=i+2,nstart_sup+nsup-1
1511           nhpb=nhpb+1
1512           ihpb(nhpb)=i+nstart_seq-nstart_sup
1513           jhpb(nhpb)=j+nstart_seq-nstart_sup
1514           forcon(nhpb)=weidis
1515           dhpb(nhpb)=dist(i,j)
1516         enddo
1517       enddo
1518 cd      write (iout,'(a)') 'Distance constraints:'
1519 cd      do i=nss+1,nhpb
1520 cd        ii=ihpb(i)
1521 cd        jj=jhpb(i)
1522 cd        iden='CA'
1523 cd        if (ii.gt.nres) then
1524 cd          iden='SC'
1525 cd          ii=ii-nres
1526 cd          jj=jj-nres
1527 cd        endif
1528 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1529 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1530 cd     &  dhpb(i),forcon(i)
1531 cd      enddo
1532       return
1533       end
1534 c----------------------------------------------------------------------------
1535       subroutine map_read
1536       implicit real*8 (a-h,o-z)
1537       include 'DIMENSIONS'
1538       include 'COMMON.MAP'
1539       include 'COMMON.IOUNITS'
1540       character*3 angid(4) /'THE','PHI','ALP','OME'/
1541       character*80 mapcard,ucase
1542       do imap=1,nmap
1543         read (inp,'(a)') mapcard
1544         mapcard=ucase(mapcard)
1545         if (index(mapcard,'PHI').gt.0) then
1546           kang(imap)=1
1547         else if (index(mapcard,'THE').gt.0) then
1548           kang(imap)=2
1549         else if (index(mapcard,'ALP').gt.0) then
1550           kang(imap)=3
1551         else if (index(mapcard,'OME').gt.0) then
1552           kang(imap)=4
1553         else
1554           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1555           stop 'Error - illegal variable spec in MAP card.'
1556         endif
1557         call readi (mapcard,'RES1',res1(imap),0)
1558         call readi (mapcard,'RES2',res2(imap),0)
1559         if (res1(imap).eq.0) then
1560           res1(imap)=res2(imap)
1561         else if (res2(imap).eq.0) then
1562           res2(imap)=res1(imap)
1563         endif
1564         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1565           write (iout,'(a)')
1566      &    'Error - illegal definition of variable group in MAP.'
1567           stop 'Error - illegal definition of variable group in MAP.'
1568         endif
1569         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1570         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1571         call readi(mapcard,'NSTEP',nstep(imap),0)
1572         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1573           write (iout,'(a)')
1574      &     'Illegal boundary and/or step size specification in MAP.'
1575           stop 'Illegal boundary and/or step size specification in MAP.'
1576         endif
1577       enddo ! imap
1578       return
1579       end
1580 c----------------------------------------------------------------------------
1581       subroutine csaread
1582       implicit real*8 (a-h,o-z)
1583       include 'DIMENSIONS'
1584       include 'COMMON.IOUNITS'
1585       include 'COMMON.GEO'
1586       include 'COMMON.CSA'
1587       include 'COMMON.BANK'
1588       include 'COMMON.CONTROL'
1589       character*80 ucase
1590       character*620 mcmcard
1591       call card_concat(mcmcard)
1592
1593       call readi(mcmcard,'NCONF',nconf,50)
1594       call readi(mcmcard,'NADD',nadd,0)
1595       call readi(mcmcard,'JSTART',jstart,1)
1596       call readi(mcmcard,'JEND',jend,1)
1597       call readi(mcmcard,'NSTMAX',nstmax,500000)
1598       call readi(mcmcard,'N0',n0,1)
1599       call readi(mcmcard,'N1',n1,6)
1600       call readi(mcmcard,'N2',n2,4)
1601       call readi(mcmcard,'N3',n3,0)
1602       call readi(mcmcard,'N4',n4,0)
1603       call readi(mcmcard,'N5',n5,0)
1604       call readi(mcmcard,'N6',n6,10)
1605       call readi(mcmcard,'N7',n7,0)
1606       call readi(mcmcard,'N8',n8,0)
1607       call readi(mcmcard,'N9',n9,0)
1608       call readi(mcmcard,'N14',n14,0)
1609       call readi(mcmcard,'N15',n15,0)
1610       call readi(mcmcard,'N16',n16,0)
1611       call readi(mcmcard,'N17',n17,0)
1612       call readi(mcmcard,'N18',n18,0)
1613
1614       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1615
1616       call readi(mcmcard,'NDIFF',ndiff,2)
1617       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1618       call readi(mcmcard,'IS1',is1,1)
1619       call readi(mcmcard,'IS2',is2,8)
1620       call readi(mcmcard,'NRAN0',nran0,4)
1621       call readi(mcmcard,'NRAN1',nran1,2)
1622       call readi(mcmcard,'IRR',irr,1)
1623       call readi(mcmcard,'NSEED',nseed,20)
1624       call readi(mcmcard,'NTOTAL',ntotal,10000)
1625       call reada(mcmcard,'CUT1',cut1,2.0d0)
1626       call reada(mcmcard,'CUT2',cut2,5.0d0)
1627       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1628       call readi(mcmcard,'ICMAX',icmax,3)
1629       call readi(mcmcard,'IRESTART',irestart,0)
1630 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1631       ntbankm=0
1632 c!bankt
1633       call reada(mcmcard,'DELE',dele,20.0d0)
1634       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1635       call readi(mcmcard,'IREF',iref,0)
1636       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1637       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1638       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1639       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1640       write (iout,*) "NCONF_IN",nconf_in
1641       return
1642       end
1643 c----------------------------------------------------------------------------
1644 cfmc      subroutine mcmfread
1645 cfmc      implicit real*8 (a-h,o-z)
1646 cfmc      include 'DIMENSIONS'
1647 cfmc      include 'COMMON.MCMF'
1648 cfmc      include 'COMMON.IOUNITS'
1649 cfmc      include 'COMMON.GEO'
1650 cfmc      character*80 ucase
1651 cfmc      character*620 mcmcard
1652 cfmc      call card_concat(mcmcard)
1653 cfmc
1654 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1655 cfmc      write(iout,*)'MAXRANT=',maxrant
1656 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1657 cfmc      write(iout,*)'MAXFAM=',maxfam
1658 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1659 cfmc      write(iout,*)'NNET1=',nnet1
1660 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1661 cfmc      write(iout,*)'NNET2=',nnet2
1662 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1663 cfmc      write(iout,*)'NNET3=',nnet3
1664 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1665 cfmc      write(iout,*)'ILASTT=',ilastt
1666 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1667 cfmc      write(iout,*)'MAXSTR=',maxstr
1668 cfmc      maxstr_f=maxstr/maxfam
1669 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1670 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1671 cfmc      write(iout,*)'NMCMF=',nmcmf
1672 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1673 cfmc      write(iout,*)'IFOCUS=',ifocus
1674 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1675 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1676 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1677 cfmc      write(iout,*)'INTPRT=',intprt
1678 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1679 cfmc      write(iout,*)'IPRT=',iprt
1680 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1681 cfmc      write(iout,*)'IMAXTR=',imaxtr
1682 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1683 cfmc      write(iout,*)'MAXEVEN=',maxeven
1684 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1685 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1686 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1687 cfmc      write(iout,*)'INIMIN=',inimin
1688 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1689 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1690 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1691 cfmc      write(iout,*)'NTHREAD=',nthread
1692 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1693 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1694 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1695 cfmc      write(iout,*)'MAXPERT=',maxpert
1696 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1697 cfmc      write(iout,*)'IRMSD=',irmsd
1698 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1699 cfmc      write(iout,*)'DENEMIN=',denemin
1700 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1701 cfmc      write(iout,*)'RCUT1S=',rcut1s
1702 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1703 cfmc      write(iout,*)'RCUT1E=',rcut1e
1704 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1705 cfmc      write(iout,*)'RCUT2S=',rcut2s
1706 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1707 cfmc      write(iout,*)'RCUT2E=',rcut2e
1708 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1709 cfmc      write(iout,*)'DPERT1=',d_pert1
1710 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1711 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1712 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1713 cfmc      write(iout,*)'DPERT2=',d_pert2
1714 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1715 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1716 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1717 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1718 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1719 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1720 cfmc      d_pert1=deg2rad*d_pert1
1721 cfmc      d_pert1a=deg2rad*d_pert1a
1722 cfmc      d_pert2=deg2rad*d_pert2
1723 cfmc      d_pert2a=deg2rad*d_pert2a
1724 cfmc      d_pert2b=deg2rad*d_pert2b
1725 cfmc      d_pert2c=deg2rad*d_pert2c
1726 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1727 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1728 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1729 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1730 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1731 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1732 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1733 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1734 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1735 cfmc      write(iout,*)'RCUTINI=',rcutini
1736 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1737 cfmc      write(iout,*)'GRAT=',grat
1738 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1739 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1740 cfmc
1741 cfmc      return
1742 cfmc      end
1743 c----------------------------------------------------------------------------
1744       subroutine mcmread
1745       implicit real*8 (a-h,o-z)
1746       include 'DIMENSIONS'
1747       include 'COMMON.MCM'
1748       include 'COMMON.MCE'
1749       include 'COMMON.IOUNITS'
1750       character*80 ucase
1751       character*320 mcmcard
1752       call card_concat(mcmcard)
1753       call readi(mcmcard,'MAXACC',maxacc,100)
1754       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1755       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1756       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1757       call readi(mcmcard,'MAXREPM',maxrepm,200)
1758       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1759       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1760       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1761       call reada(mcmcard,'E_UP',e_up,5.0D0)
1762       call reada(mcmcard,'DELTE',delte,0.1D0)
1763       call readi(mcmcard,'NSWEEP',nsweep,5)
1764       call readi(mcmcard,'NSTEPH',nsteph,0)
1765       call readi(mcmcard,'NSTEPC',nstepc,0)
1766       call reada(mcmcard,'TMIN',tmin,298.0D0)
1767       call reada(mcmcard,'TMAX',tmax,298.0D0)
1768       call readi(mcmcard,'NWINDOW',nwindow,0)
1769       call readi(mcmcard,'PRINT_MC',print_mc,0)
1770       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1771       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1772       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1773       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1774       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1775       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1776       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1777       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1778       if (nwindow.gt.0) then
1779         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1780         do i=1,nwindow
1781           winlen(i)=winend(i)-winstart(i)+1
1782         enddo
1783       endif
1784       if (tmax.lt.tmin) tmax=tmin
1785       if (tmax.eq.tmin) then
1786         nstepc=0
1787         nsteph=0
1788       endif
1789       if (nstepc.gt.0 .and. nsteph.gt.0) then
1790         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1791         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1792       endif
1793 C Probabilities of different move types
1794       sumpro_type(0)=0.0D0
1795       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1796       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1797       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1798       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1799       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1800       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1801       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1802       do i=1,MaxMoveType
1803         print *,'i',i,' sumprotype',sumpro_type(i)
1804         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1805         print *,'i',i,' sumprotype',sumpro_type(i)
1806       enddo
1807       return
1808       end
1809 c----------------------------------------------------------------------------
1810       subroutine read_minim
1811       implicit real*8 (a-h,o-z)
1812       include 'DIMENSIONS'
1813       include 'COMMON.MINIM'
1814       include 'COMMON.IOUNITS'
1815       character*80 ucase
1816       character*320 minimcard
1817       call card_concat(minimcard)
1818       call readi(minimcard,'MAXMIN',maxmin,2000)
1819       call readi(minimcard,'MAXFUN',maxfun,5000)
1820       call readi(minimcard,'MINMIN',minmin,maxmin)
1821       call readi(minimcard,'MINFUN',minfun,maxmin)
1822       call reada(minimcard,'TOLF',tolf,1.0D-2)
1823       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1824       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1825       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1826       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1827       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1828      &         'Options in energy minimization:'
1829       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1830      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1831      & 'MinMin:',MinMin,' MinFun:',MinFun,
1832      & ' TolF:',TolF,' RTolF:',RTolF
1833       return
1834       end
1835 c----------------------------------------------------------------------------
1836       subroutine read_angles(kanal,*)
1837       implicit real*8 (a-h,o-z)
1838       include 'DIMENSIONS'
1839       include 'COMMON.GEO'
1840       include 'COMMON.VAR'
1841       include 'COMMON.CHAIN'
1842       include 'COMMON.IOUNITS'
1843       include 'COMMON.CONTROL'
1844 c Read angles from input
1845 c
1846        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1847        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1848        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1849        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1850
1851        do i=1,nres
1852 c 9/7/01 avoid 180 deg valence angle
1853         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1854 c
1855         theta(i)=deg2rad*theta(i)
1856         phi(i)=deg2rad*phi(i)
1857         alph(i)=deg2rad*alph(i)
1858         omeg(i)=deg2rad*omeg(i)
1859        enddo
1860       return
1861    10 return1
1862       end
1863 c----------------------------------------------------------------------------
1864       subroutine reada(rekord,lancuch,wartosc,default)
1865       implicit none
1866       character*(*) rekord,lancuch
1867       double precision wartosc,default
1868       integer ilen,iread
1869       external ilen
1870       iread=index(rekord,lancuch)
1871       if (iread.eq.0) then
1872         wartosc=default
1873         return
1874       endif
1875       iread=iread+ilen(lancuch)+1
1876       read (rekord(iread:),*,err=10,end=10) wartosc
1877       return
1878   10  wartosc=default
1879       return
1880       end
1881 c----------------------------------------------------------------------------
1882       subroutine readi(rekord,lancuch,wartosc,default)
1883       implicit none
1884       character*(*) rekord,lancuch
1885       integer wartosc,default
1886       integer ilen,iread
1887       external ilen
1888       iread=index(rekord,lancuch)
1889       if (iread.eq.0) then
1890         wartosc=default
1891         return
1892       endif
1893       iread=iread+ilen(lancuch)+1
1894       read (rekord(iread:),*,err=10,end=10) wartosc
1895       return
1896   10  wartosc=default
1897       return
1898       end
1899 c----------------------------------------------------------------------------
1900       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1901       implicit none
1902       integer dim,i
1903       integer tablica(dim),default
1904       character*(*) rekord,lancuch
1905       character*80 aux
1906       integer ilen,iread
1907       external ilen
1908       do i=1,dim
1909         tablica(i)=default
1910       enddo
1911       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1912       if (iread.eq.0) return
1913       iread=iread+ilen(lancuch)+1
1914       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1915    10 return
1916       end
1917 c----------------------------------------------------------------------------
1918       subroutine multreada(rekord,lancuch,tablica,dim,default)
1919       implicit none
1920       integer dim,i
1921       double precision tablica(dim),default
1922       character*(*) rekord,lancuch
1923       character*80 aux
1924       integer ilen,iread
1925       external ilen
1926       do i=1,dim
1927         tablica(i)=default
1928       enddo
1929       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1930       if (iread.eq.0) return
1931       iread=iread+ilen(lancuch)+1
1932       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1933    10 return
1934       end
1935 c----------------------------------------------------------------------------
1936       subroutine openunits
1937       implicit real*8 (a-h,o-z)
1938       include 'DIMENSIONS'
1939 #ifdef MPI
1940       include 'mpif.h'
1941       character*16 form,nodename
1942       integer nodelen
1943 #endif
1944       include 'COMMON.SETUP'
1945       include 'COMMON.IOUNITS'
1946       include 'COMMON.MD'
1947       include 'COMMON.CONTROL'
1948       integer lenpre,lenpot,ilen,lentmp
1949       external ilen
1950       character*3 out1file_text,ucase
1951       character*3 ll
1952       external ucase
1953 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1954       call getenv_loc("PREFIX",prefix)
1955       pref_orig = prefix
1956       call getenv_loc("POT",pot)
1957       call getenv_loc("DIRTMP",tmpdir)
1958       call getenv_loc("CURDIR",curdir)
1959       call getenv_loc("OUT1FILE",out1file_text)
1960 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1961       out1file_text=ucase(out1file_text)
1962       if (out1file_text(1:1).eq."Y") then
1963         out1file=.true.
1964       else
1965         out1file=fg_rank.gt.0
1966       endif
1967       lenpre=ilen(prefix)
1968       lenpot=ilen(pot)
1969       lentmp=ilen(tmpdir)
1970       if (lentmp.gt.0) then
1971           write (*,'(80(1h!))')
1972           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1973           write (*,'(80(1h!))')
1974           write (*,*)"All output files will be on node /tmp directory."
1975 #ifdef MPI
1976         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1977         if (me.eq.king) then
1978           write (*,*) "The master node is ",nodename
1979         else if (fg_rank.eq.0) then
1980           write (*,*) "I am the CG slave node ",nodename
1981         else
1982           write (*,*) "I am the FG slave node ",nodename
1983         endif
1984 #endif
1985         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1986         lenpre = lentmp+lenpre+1
1987       endif
1988       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1989 C Get the names and open the input files
1990 #if defined(WINIFL) || defined(WINPGI)
1991       open(1,file=pref_orig(:ilen(pref_orig))//
1992      &  '.inp',status='old',readonly,shared)
1993        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1994 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1995 C Get parameter filenames and open the parameter files.
1996       call getenv_loc('BONDPAR',bondname)
1997       open (ibond,file=bondname,status='old',readonly,shared)
1998       call getenv_loc('THETPAR',thetname)
1999       open (ithep,file=thetname,status='old',readonly,shared)
2000       call getenv_loc('ROTPAR',rotname)
2001       open (irotam,file=rotname,status='old',readonly,shared)
2002       call getenv_loc('TORPAR',torname)
2003       open (itorp,file=torname,status='old',readonly,shared)
2004       call getenv_loc('TORDPAR',tordname)
2005       open (itordp,file=tordname,status='old',readonly,shared)
2006       call getenv_loc('TORKCC',torkccname)
2007       open (itorkcc,file=torkccname,status='old',readonly,shared)
2008       call getenv_loc('THETKCC',thetkccname)
2009       open (ithetkcc,file=thetkccname,status='old',readonly,shared)
2010       call getenv_loc('FOURIER',fouriername)
2011       open (ifourier,file=fouriername,status='old',readonly,shared)
2012       call getenv_loc('ELEPAR',elename)
2013       open (ielep,file=elename,status='old',readonly,shared)
2014       call getenv_loc('SIDEPAR',sidename)
2015       open (isidep,file=sidename,status='old',readonly,shared)
2016       call getenv_loc('LIPTRANPAR',liptranname)
2017       open (iliptranpar,file=liptranname,status='old',readonly,shared)
2018 #elif (defined CRAY) || (defined AIX)
2019       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2020      &  action='read')
2021 c      print *,"Processor",myrank," opened file 1"
2022       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2023 c      print *,"Processor",myrank," opened file 9"
2024 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2025 C Get parameter filenames and open the parameter files.
2026       call getenv_loc('BONDPAR',bondname)
2027       open (ibond,file=bondname,status='old',action='read')
2028 c      print *,"Processor",myrank," opened file IBOND"
2029       call getenv_loc('THETPAR',thetname)
2030       open (ithep,file=thetname,status='old',action='read')
2031 c      print *,"Processor",myrank," opened file ITHEP"
2032       call getenv_loc('ROTPAR',rotname)
2033       open (irotam,file=rotname,status='old',action='read')
2034 c      print *,"Processor",myrank," opened file IROTAM"
2035       call getenv_loc('TORPAR',torname)
2036       open (itorp,file=torname,status='old',action='read')
2037 c      print *,"Processor",myrank," opened file ITORP"
2038       call getenv_loc('TORDPAR',tordname)
2039       open (itordp,file=tordname,status='old',action='read')
2040       call getenv_loc('TORKCC',torkccname)
2041       open (itorkcc,file=torkccname,status='old',action='read')
2042       call getenv_loc('THETKCC',thetkccname)
2043       open (ithetkcc,file=thetkccname,status='old',action='read')
2044 c      print *,"Processor",myrank," opened file ITORDP"
2045       call getenv_loc('SCCORPAR',sccorname)
2046       open (isccor,file=sccorname,status='old',action='read')
2047 c      print *,"Processor",myrank," opened file ISCCOR"
2048       call getenv_loc('FOURIER',fouriername)
2049       open (ifourier,file=fouriername,status='old',action='read')
2050 c      print *,"Processor",myrank," opened file IFOURIER"
2051       call getenv_loc('ELEPAR',elename)
2052       open (ielep,file=elename,status='old',action='read')
2053 c      print *,"Processor",myrank," opened file IELEP"
2054       call getenv_loc('SIDEPAR',sidename)
2055       open (isidep,file=sidename,status='old',action='read')
2056       call getenv_loc('LIPTRANPAR',liptranname)
2057       open (iliptranpar,file=liptranname,status='old',action='read')
2058 c      print *,"Processor",myrank," opened file ISIDEP"
2059 c      print *,"Processor",myrank," opened parameter files"
2060 #elif (defined G77)
2061       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2062       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2063 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2064 C Get parameter filenames and open the parameter files.
2065       call getenv_loc('BONDPAR',bondname)
2066       open (ibond,file=bondname,status='old')
2067       call getenv_loc('THETPAR',thetname)
2068       open (ithep,file=thetname,status='old')
2069       call getenv_loc('ROTPAR',rotname)
2070       open (irotam,file=rotname,status='old')
2071       call getenv_loc('TORPAR',torname)
2072       open (itorp,file=torname,status='old')
2073       call getenv_loc('TORDPAR',tordname)
2074       open (itordp,file=tordname,status='old')
2075       call getenv_loc('TORKCC',torkccname)
2076       open (itorkcc,file=torkccname,status='old')
2077       call getenv_loc('THETKCC',thetkccname)
2078       open (ithetkcc,file=thetkccname,status='old')
2079       call getenv_loc('SCCORPAR',sccorname)
2080       open (isccor,file=sccorname,status='old')
2081       call getenv_loc('FOURIER',fouriername)
2082       open (ifourier,file=fouriername,status='old')
2083       call getenv_loc('ELEPAR',elename)
2084       open (ielep,file=elename,status='old')
2085       call getenv_loc('SIDEPAR',sidename)
2086       open (isidep,file=sidename,status='old')
2087       call getenv_loc('LIPTRANPAR',liptranname)
2088       open (iliptranpar,file=liptranname,status='old')
2089 #else
2090       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2091      &  readonly)
2092        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2093 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2094 C Get parameter filenames and open the parameter files.
2095       call getenv_loc('BONDPAR',bondname)
2096       open (ibond,file=bondname,status='old',readonly)
2097       call getenv_loc('THETPAR',thetname)
2098       open (ithep,file=thetname,status='old',readonly)
2099       call getenv_loc('ROTPAR',rotname)
2100       open (irotam,file=rotname,status='old',readonly)
2101       call getenv_loc('TORPAR',torname)
2102       open (itorp,file=torname,status='old',readonly)
2103       call getenv_loc('TORDPAR',tordname)
2104       open (itordp,file=tordname,status='old',readonly)
2105       call getenv_loc('TORKCC',torkccname)
2106       open (itorkcc,file=torkccname,status='old',readonly)
2107       call getenv_loc('THETKCC',thetkccname)
2108       open (ithetkcc,file=thetkccname,status='old',readonly)
2109       call getenv_loc('SCCORPAR',sccorname)
2110       open (isccor,file=sccorname,status='old',readonly)
2111 #ifndef CRYST_THETA
2112       call getenv_loc('THETPARPDB',thetname_pdb)
2113       print *,"thetname_pdb ",thetname_pdb
2114       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2115       print *,ithep_pdb," opened"
2116 #endif
2117       call getenv_loc('FOURIER',fouriername)
2118       open (ifourier,file=fouriername,status='old',readonly)
2119       call getenv_loc('ELEPAR',elename)
2120       open (ielep,file=elename,status='old',readonly)
2121       call getenv_loc('SIDEPAR',sidename)
2122       open (isidep,file=sidename,status='old',readonly)
2123       call getenv_loc('LIPTRANPAR',liptranname)
2124       open (iliptranpar,file=liptranname,status='old',action='read')
2125 #ifndef CRYST_SC
2126       call getenv_loc('ROTPARPDB',rotname_pdb)
2127       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2128 #endif
2129 #endif
2130 #ifndef OLDSCP
2131 C
2132 C 8/9/01 In the newest version SCp interaction constants are read from a file
2133 C Use -DOLDSCP to use hard-coded constants instead.
2134 C
2135       call getenv_loc('SCPPAR',scpname)
2136 #if defined(WINIFL) || defined(WINPGI)
2137       open (iscpp,file=scpname,status='old',readonly,shared)
2138 #elif (defined CRAY)  || (defined AIX)
2139       open (iscpp,file=scpname,status='old',action='read')
2140 #elif (defined G77)
2141       open (iscpp,file=scpname,status='old')
2142 #else
2143       open (iscpp,file=scpname,status='old',readonly)
2144 #endif
2145 #endif
2146       call getenv_loc('PATTERN',patname)
2147 #if defined(WINIFL) || defined(WINPGI)
2148       open (icbase,file=patname,status='old',readonly,shared)
2149 #elif (defined CRAY)  || (defined AIX)
2150       open (icbase,file=patname,status='old',action='read')
2151 #elif (defined G77)
2152       open (icbase,file=patname,status='old')
2153 #else
2154       open (icbase,file=patname,status='old',readonly)
2155 #endif
2156 #ifdef MPI
2157 C Open output file only for CG processes
2158 c      print *,"Processor",myrank," fg_rank",fg_rank
2159       if (fg_rank.eq.0) then
2160
2161       if (nodes.eq.1) then
2162         npos=3
2163       else
2164         npos = dlog10(dfloat(nodes-1))+1
2165       endif
2166       if (npos.lt.3) npos=3
2167       write (liczba,'(i1)') npos
2168       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2169      &  //')'
2170       write (liczba,form) me
2171       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2172      &  liczba(:ilen(liczba))
2173       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2174      &  //'.int'
2175       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2176      &  //'.pdb'
2177       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2178      &  liczba(:ilen(liczba))//'.mol2'
2179       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2180      &  liczba(:ilen(liczba))//'.stat'
2181       if (lentmp.gt.0)
2182      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2183      &      //liczba(:ilen(liczba))//'.stat')
2184       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2185      &  //'.rst'
2186       if(usampl) then
2187           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2188      & liczba(:ilen(liczba))//'.const'
2189       endif
2190
2191       endif
2192 #else
2193       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2194       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2195       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2196       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2197       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2198       if (lentmp.gt.0)
2199      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2200      &    //'.stat')
2201       rest2name=prefix(:ilen(prefix))//'.rst'
2202       if(usampl) then
2203          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2204       endif
2205 #endif
2206 #if defined(AIX) || defined(PGI)
2207       if (me.eq.king .or. .not. out1file)
2208      &   open(iout,file=outname,status='unknown')
2209 #ifdef DEBUG
2210       if (fg_rank.gt.0) then
2211         write (liczba,'(i3.3)') myrank/nfgtasks
2212         write (ll,'(bz,i3.3)') fg_rank
2213         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2214      &   status='unknown')
2215       endif
2216 #endif
2217       if(me.eq.king) then
2218        open(igeom,file=intname,status='unknown',position='append')
2219        open(ipdb,file=pdbname,status='unknown')
2220        open(imol2,file=mol2name,status='unknown')
2221        open(istat,file=statname,status='unknown',position='append')
2222       else
2223 c1out       open(iout,file=outname,status='unknown')
2224       endif
2225 #else
2226       if (me.eq.king .or. .not.out1file)
2227      &    open(iout,file=outname,status='unknown')
2228 #ifdef DEBUG
2229       if (fg_rank.gt.0) then
2230         write (liczba,'(i3.3)') myrank/nfgtasks
2231         write (ll,'(bz,i3.3)') fg_rank
2232         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2233      &   status='unknown')
2234       endif
2235 #endif
2236       if(me.eq.king) then
2237        open(igeom,file=intname,status='unknown',access='append')
2238        open(ipdb,file=pdbname,status='unknown')
2239        open(imol2,file=mol2name,status='unknown')
2240        open(istat,file=statname,status='unknown',access='append')
2241       else
2242 c1out       open(iout,file=outname,status='unknown')
2243       endif
2244 #endif
2245       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2246       csa_seed=prefix(:lenpre)//'.CSA.seed'
2247       csa_history=prefix(:lenpre)//'.CSA.history'
2248       csa_bank=prefix(:lenpre)//'.CSA.bank'
2249       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2250       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2251       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2252 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2253       csa_int=prefix(:lenpre)//'.int'
2254       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2255       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2256       csa_in=prefix(:lenpre)//'.CSA.in'
2257 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2258 C Write file names
2259       if (me.eq.king)then
2260       write (iout,'(80(1h-))')
2261       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2262       write (iout,'(80(1h-))')
2263       write (iout,*) "Input file                      : ",
2264      &  pref_orig(:ilen(pref_orig))//'.inp'
2265       write (iout,*) "Output file                     : ",
2266      &  outname(:ilen(outname))
2267       write (iout,*)
2268       write (iout,*) "Sidechain potential file        : ",
2269      &  sidename(:ilen(sidename))
2270 #ifndef OLDSCP
2271       write (iout,*) "SCp potential file              : ",
2272      &  scpname(:ilen(scpname))
2273 #endif
2274       write (iout,*) "Electrostatic potential file    : ",
2275      &  elename(:ilen(elename))
2276       write (iout,*) "Cumulant coefficient file       : ",
2277      &  fouriername(:ilen(fouriername))
2278       write (iout,*) "Torsional parameter file        : ",
2279      &  torname(:ilen(torname))
2280       write (iout,*) "Double torsional parameter file : ",
2281      &  tordname(:ilen(tordname))
2282       write (iout,*) "SCCOR parameter file : ",
2283      &  sccorname(:ilen(sccorname))
2284       write (iout,*) "Bond & inertia constant file    : ",
2285      &  bondname(:ilen(bondname))
2286       write (iout,*) "Bending parameter file          : ",
2287      &  thetname(:ilen(thetname))
2288       write (iout,*) "Rotamer parameter file          : ",
2289      &  rotname(:ilen(rotname))
2290       write (iout,*) "Thetpdb parameter file          : ",
2291      &  thetname_pdb(:ilen(thetname_pdb))
2292       write (iout,*) "Threading database              : ",
2293      &  patname(:ilen(patname))
2294       if (lentmp.ne.0)
2295      &write (iout,*)" DIRTMP                          : ",
2296      &  tmpdir(:lentmp)
2297       write (iout,'(80(1h-))')
2298       endif
2299       return
2300       end
2301 c----------------------------------------------------------------------------
2302       subroutine card_concat(card)
2303       implicit real*8 (a-h,o-z)
2304       include 'DIMENSIONS'
2305       include 'COMMON.IOUNITS'
2306       character*(*) card
2307       character*80 karta,ucase
2308       external ilen
2309       read (inp,'(a)') karta
2310       karta=ucase(karta)
2311       card=' '
2312       do while (karta(80:80).eq.'&')
2313         card=card(:ilen(card)+1)//karta(:79)
2314         read (inp,'(a)') karta
2315         karta=ucase(karta)
2316       enddo
2317       card=card(:ilen(card)+1)//karta
2318       return
2319       end
2320 c----------------------------------------------------------------------------------
2321       subroutine readrst
2322       implicit real*8 (a-h,o-z)
2323       include 'DIMENSIONS'
2324       include 'COMMON.CHAIN'
2325       include 'COMMON.IOUNITS'
2326       include 'COMMON.MD'
2327       open(irest2,file=rest2name,status='unknown')
2328       read(irest2,*) totT,EK,potE,totE,t_bath
2329       totTafm=totT
2330       do i=1,2*nres
2331          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2332       enddo
2333       do i=1,2*nres
2334          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2335       enddo
2336       if(usampl) then
2337              read (irest2,*) iset
2338       endif
2339       close(irest2)
2340       return
2341       end
2342 c---------------------------------------------------------------------------------
2343       subroutine read_fragments
2344       implicit real*8 (a-h,o-z)
2345       include 'DIMENSIONS'
2346 #ifdef MPI
2347       include 'mpif.h'
2348 #endif
2349       include 'COMMON.SETUP'
2350       include 'COMMON.CHAIN'
2351       include 'COMMON.IOUNITS'
2352       include 'COMMON.MD'
2353       include 'COMMON.CONTROL'
2354       read(inp,*) nset,nfrag,npair,nfrag_back
2355       if(me.eq.king.or..not.out1file)
2356      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2357      &  " nfrag_back",nfrag_back
2358       do iset=1,nset
2359          read(inp,*) mset(iset)
2360        do i=1,nfrag
2361          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2362      &     qinfrag(i,iset)
2363          if(me.eq.king.or..not.out1file)
2364      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2365      &     ifrag(2,i,iset), qinfrag(i,iset)
2366        enddo
2367        do i=1,npair
2368         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2369      &    qinpair(i,iset)
2370         if(me.eq.king.or..not.out1file)
2371      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2372      &    ipair(2,i,iset), qinpair(i,iset)
2373        enddo
2374        do i=1,nfrag_back
2375         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2376      &     wfrag_back(3,i,iset),
2377      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2378         if(me.eq.king.or..not.out1file)
2379      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2380      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2381        enddo
2382       enddo
2383       return
2384       end
2385 C---------------------------------------------------------------------------
2386       subroutine read_afminp
2387             implicit real*8 (a-h,o-z)
2388       include 'DIMENSIONS'
2389 #ifdef MPI
2390       include 'mpif.h'
2391 #endif
2392       include 'COMMON.SETUP'
2393       include 'COMMON.CONTROL'
2394       include 'COMMON.CHAIN'
2395       include 'COMMON.IOUNITS'
2396       include 'COMMON.SBRIDGE'
2397       character*320 afmcard
2398       print *, "wchodze"
2399       call card_concat(afmcard)
2400       call readi(afmcard,"BEG",afmbeg,0)
2401       call readi(afmcard,"END",afmend,0)
2402       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2403       call reada(afmcard,"VEL",velAFMconst,0.0d0)
2404       print *,'FORCE=' ,forceAFMconst
2405 CCCC NOW PROPERTIES FOR AFM
2406        distafminit=0.0d0
2407        do i=1,3
2408         distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2409        enddo
2410         distafminit=dsqrt(distafminit)
2411         print *,'initdist',distafminit
2412       return
2413       end
2414
2415 c-------------------------------------------------------------------------------
2416       subroutine read_dist_constr
2417       implicit real*8 (a-h,o-z)
2418       include 'DIMENSIONS'
2419 #ifdef MPI
2420       include 'mpif.h'
2421 #endif
2422       include 'COMMON.SETUP'
2423       include 'COMMON.CONTROL'
2424       include 'COMMON.CHAIN'
2425       include 'COMMON.IOUNITS'
2426       include 'COMMON.SBRIDGE'
2427       integer ifrag_(2,100),ipair_(2,100)
2428       double precision wfrag_(100),wpair_(100)
2429       character*500 controlcard
2430       print *, "WCHODZE"
2431       write (iout,*) "Calling read_dist_constr"
2432 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2433 c      call flush(iout)
2434       call card_concat(controlcard)
2435       call readi(controlcard,"NFRAG",nfrag_,0)
2436       call readi(controlcard,"NPAIR",npair_,0)
2437       call readi(controlcard,"NDIST",ndist_,0)
2438       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2439       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2440       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2441       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2442       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2443 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2444 c      write (iout,*) "IFRAG"
2445 c      do i=1,nfrag_
2446 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2447 c      enddo
2448 c      write (iout,*) "IPAIR"
2449 c      do i=1,npair_
2450 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2451 c      enddo
2452       call flush(iout)
2453       do i=1,nfrag_
2454         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2455         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2456      &    ifrag_(2,i)=nstart_sup+nsup-1
2457 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2458         call flush(iout)
2459         if (wfrag_(i).gt.0.0d0) then
2460         do j=ifrag_(1,i),ifrag_(2,i)-1
2461           do k=j+1,ifrag_(2,i)
2462 c            write (iout,*) "j",j," k",k
2463             ddjk=dist(j,k)
2464             if (constr_dist.eq.1) then
2465             nhpb=nhpb+1
2466             ihpb(nhpb)=j
2467             jhpb(nhpb)=k
2468               dhpb(nhpb)=ddjk
2469             forcon(nhpb)=wfrag_(i)
2470             else if (constr_dist.eq.2) then
2471               if (ddjk.le.dist_cut) then
2472                 nhpb=nhpb+1
2473                 ihpb(nhpb)=j
2474                 jhpb(nhpb)=k
2475                 dhpb(nhpb)=ddjk
2476                 forcon(nhpb)=wfrag_(i)
2477               endif
2478             else
2479               nhpb=nhpb+1
2480               ihpb(nhpb)=j
2481               jhpb(nhpb)=k
2482               dhpb(nhpb)=ddjk
2483               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2484             endif
2485 #ifdef MPI
2486             if (.not.out1file .or. me.eq.king)
2487      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2488      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2489 #else
2490             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2491      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2492 #endif
2493           enddo
2494         enddo
2495         endif
2496       enddo
2497       do i=1,npair_
2498         if (wpair_(i).gt.0.0d0) then
2499         ii = ipair_(1,i)
2500         jj = ipair_(2,i)
2501         if (ii.gt.jj) then
2502           itemp=ii
2503           ii=jj
2504           jj=itemp
2505         endif
2506         do j=ifrag_(1,ii),ifrag_(2,ii)
2507           do k=ifrag_(1,jj),ifrag_(2,jj)
2508             nhpb=nhpb+1
2509             ihpb(nhpb)=j
2510             jhpb(nhpb)=k
2511             forcon(nhpb)=wpair_(i)
2512             dhpb(nhpb)=dist(j,k)
2513 #ifdef MPI
2514             if (.not.out1file .or. me.eq.king)
2515      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2516      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2517 #else
2518             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2519      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2520 #endif
2521           enddo
2522         enddo
2523         endif
2524       enddo
2525       print *,ndist_
2526       do i=1,ndist_
2527         if (constr_dist.eq.11) then
2528         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2529      &     ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2530         fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2531         else
2532 C        print *,"in else"
2533         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2534      &     ibecarb(i),forcon(nhpb+1)
2535         endif
2536         if (forcon(nhpb+1).gt.0.0d0) then
2537           nhpb=nhpb+1
2538           if (ibecarb(i).gt.0) then
2539             ihpb(i)=ihpb(i)+nres
2540             jhpb(i)=jhpb(i)+nres
2541           endif
2542           if (dhpb(nhpb).eq.0.0d0)
2543      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2544         endif
2545 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2546 C        if (forcon(nhpb+1).gt.0.0d0) then
2547 C          nhpb=nhpb+1
2548 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2549 #ifdef MPI
2550           if (.not.out1file .or. me.eq.king)
2551      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2552      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2553 #else
2554           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2555      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2556 #endif
2557
2558       enddo
2559       call flush(iout)
2560       return
2561       end
2562 c-------------------------------------------------------------------------------
2563 #ifdef WINIFL
2564       subroutine flush(iu)
2565       return
2566       end
2567 #endif
2568 #ifdef AIX
2569       subroutine flush(iu)
2570       call flush_(iu)
2571       return
2572       end
2573 #endif
2574 c------------------------------------------------------------------------------
2575       subroutine copy_to_tmp(source)
2576       include "DIMENSIONS"
2577       include "COMMON.IOUNITS"
2578       character*(*) source
2579       character* 256 tmpfile
2580       integer ilen
2581       external ilen
2582       logical ex
2583       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2584       inquire(file=tmpfile,exist=ex)
2585       if (ex) then
2586         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2587      &   " to temporary directory..."
2588         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2589         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2590       endif
2591       return
2592       end
2593 c------------------------------------------------------------------------------
2594       subroutine move_from_tmp(source)
2595       include "DIMENSIONS"
2596       include "COMMON.IOUNITS"
2597       character*(*) source
2598       integer ilen
2599       external ilen
2600       write (*,*) "Moving ",source(:ilen(source)),
2601      & " from temporary directory to working directory"
2602       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2603       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2604       return
2605       end
2606 c------------------------------------------------------------------------------
2607       subroutine random_init(seed)
2608 C
2609 C Initialize random number generator
2610 C
2611       implicit real*8 (a-h,o-z)
2612       include 'DIMENSIONS'
2613 #ifdef MPI
2614       include 'mpif.h'
2615       logical OKRandom, prng_restart
2616       real*8  r1
2617       integer iseed_array(4)
2618 #endif
2619       include 'COMMON.IOUNITS'
2620       include 'COMMON.TIME1'
2621       include 'COMMON.THREAD'
2622       include 'COMMON.SBRIDGE'
2623       include 'COMMON.CONTROL'
2624       include 'COMMON.MCM'
2625       include 'COMMON.MAP'
2626       include 'COMMON.HEADER'
2627       include 'COMMON.CSA'
2628       include 'COMMON.CHAIN'
2629       include 'COMMON.MUCA'
2630       include 'COMMON.MD'
2631       include 'COMMON.FFIELD'
2632       include 'COMMON.SETUP'
2633       iseed=-dint(dabs(seed))
2634       if (iseed.eq.0) then
2635         write (iout,'(/80(1h*)/20x,a/80(1h*))')
2636      &    'Random seed undefined. The program will stop.'
2637         write (*,'(/80(1h*)/20x,a/80(1h*))')
2638      &    'Random seed undefined. The program will stop.'
2639 #ifdef MPI
2640         call mpi_finalize(mpi_comm_world,ierr)
2641 #endif
2642         stop 'Bad random seed.'
2643       endif
2644 #ifdef MPI
2645       if (fg_rank.eq.0) then
2646       seed=seed*(me+1)+1
2647 #ifdef AMD64
2648       if(me.eq.king)
2649      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2650       OKRandom = prng_restart(me,iseed)
2651 #else
2652       do i=1,4
2653        tmp=65536.0d0**(4-i)
2654        iseed_array(i) = dint(seed/tmp)
2655        seed=seed-iseed_array(i)*tmp
2656       enddo
2657       if(me.eq.king)
2658      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2659      &                 (iseed_array(i),i=1,4)
2660       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2661      &                 (iseed_array(i),i=1,4)
2662       OKRandom = prng_restart(me,iseed_array)
2663 #endif
2664       if (OKRandom) then
2665 c        r1 = prng_next(me)
2666         r1=ran_number(0.0D0,1.0D0)
2667         if(me.eq.king)
2668      &   write (iout,*) 'ran_num',r1
2669         if (r1.lt.0.0d0) OKRandom=.false.
2670       endif
2671       if (.not.OKRandom) then
2672         write (iout,*) 'PRNG IS NOT WORKING!!!'
2673         print *,'PRNG IS NOT WORKING!!!'
2674         if (me.eq.0) then
2675          call flush(iout)
2676          call mpi_abort(mpi_comm_world,error_msg,ierr)
2677          stop
2678         else
2679          write (iout,*) 'too many processors for parallel prng'
2680          write (*,*) 'too many processors for parallel prng'
2681          call flush(iout)
2682          stop
2683         endif
2684       endif
2685       endif
2686 #else
2687       call vrndst(iseed)
2688       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
2689 #endif
2690       return
2691       end