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