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