source/unres/src_CSA/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       logical file_exist
  12 C Read force-field parameters except weights
  13       call parmread
  14 C Read job setup parameters
  15       call read_control
  16 C Read control parameters for energy minimzation if required
  17       if (minim) call read_minim
  18 C Read MCM control parameters if required
  19 c      if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  20 C Read MD control parameters if reqjuired
  21 c      if (modecalc.eq.12) call read_MDpar
  22 C Read MREMD control parameters if required
  23 c      if (modecalc.eq.14) then
  24 c         call read_MDpar
  25 c         call read_REMDpar
  26 c      endif
  27 C Read MUCA control parameters if required
  28 c      if (lmuca) call read_muca
  29 C Read CSA control parameters if required (from fort.40 if exists
  30 C otherwise from general input file)
  31       if (modecalc.eq.8) then
  32        inquire (file="fort.40",exist=file_exist)
  33        if (.not.file_exist) call csaread
  34       endif
  35 cfmc      if (modecalc.eq.10) call mcmfread
  36 C Read molecule information, molecule geometry, energy-term weights, and
  37 C restraints if requested
  38       call molread
  39
  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 c      include 'COMMON.THREAD'
  71       include 'COMMON.SBRIDGE'
  72       include 'COMMON.CONTROL'
  73       include 'COMMON.MCM'
  74 c      include 'COMMON.MAP'
  75       include 'COMMON.HEADER'
  76       include 'COMMON.CSA'
  77       include 'COMMON.CHAIN'
  78 c      include 'COMMON.MUCA'
  79 c      include 'COMMON.MD'
  80       include 'COMMON.FFIELD'
  81       include 'COMMON.SETUP'
  82       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
  83       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
  84       character*80 ucase
  85       character*320 controlcard
  86
  87       nglob_csa=0
  88       eglob_csa=1d99
  89       nmin_csa=0
  90       read (INP,'(a)') titel
  91       call card_concat(controlcard)
  92 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
  93 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
  94       call reada(controlcard,'SEED',seed,0.0D0)
  95       call random_init(seed)
  96 C Set up the time limit (caution! The time must be input in minutes!)
  97       read_cart=index(controlcard,'READ_CART').gt.0
  98       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  99       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
 100       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 101       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 102       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 103       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 104       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 105       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 106       call reada(controlcard,'DRMS',drms,0.1D0)
 107       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 108        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 109        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 110        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 111        write (iout,'(a,f10.1)')'DRMS    = ',drms
 112        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 113        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 114       endif
 115       call readi(controlcard,'NZ_START',nz_start,0)
 116       call readi(controlcard,'NZ_END',nz_end,0)
 117       call readi(controlcard,'IZ_SC',iz_sc,0)
 118       timlim=60.0D0*timlim
 119       safety = 60.0d0*safety
 120       timem=timlim
 121       modecalc=0
 122       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 123       minim=(index(controlcard,'MINIMIZE').gt.0)
 124       dccart=(index(controlcard,'CART').gt.0)
 125       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 126       overlapsc=.not.overlapsc
 127       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
 128       searchsc=.not.searchsc
 129       sideadd=(index(controlcard,'SIDEADD').gt.0)
 130       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 131       outpdb=(index(controlcard,'PDBOUT').gt.0)
 132       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 133       pdbref=(index(controlcard,'PDBREF').gt.0)
 134       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 135       indpdb=index(controlcard,'PDBSTART')
 136       extconf=(index(controlcard,'EXTCONF').gt.0)
 137       call readi(controlcard,'IPRINT',iprint,0)
 138       call readi(controlcard,'MAXGEN',maxgen,10000)
 139       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 140       call readi(controlcard,"KDIAG",kdiag,0)
 141       call readi(controlcard,"RESCALE_MODE",rescale_mode,0)
 142       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 143      & write (iout,*) "RESCALE_MODE",rescale_mode
 144       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 145       if (index(controlcard,'REGULAR').gt.0.0D0) then
 146         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 147         modecalc=1
 148         refstr=.true.
 149       endif
 150       if (index(controlcard,'CHECKGRAD').gt.0) then
 151         modecalc=5
 152         if (index(controlcard,'CART').gt.0) then
 153           icheckgrad=1
 154         elseif (index(controlcard,'CARINT').gt.0) then
 155           icheckgrad=2
 156         else
 157           icheckgrad=3
 158         endif
 159       elseif (index(controlcard,'THREAD').gt.0) then
 160         modecalc=2
 161         call readi(controlcard,'THREAD',nthread,0)
 162         if (nthread.gt.0) then
 163           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 164         else
 165           if (fg_rank.eq.0)
 166      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 167           stop 'Error termination in Read_Control.'
 168         endif
 169       else if (index(controlcard,'MCMA').gt.0) then
 170         modecalc=3
 171       else if (index(controlcard,'MCEE').gt.0) then
 172         modecalc=6
 173       else if (index(controlcard,'MULTCONF').gt.0) then
 174         modecalc=4
 175       else if (index(controlcard,'MAP').gt.0) then
 176         modecalc=7
 177         call readi(controlcard,'MAP',nmap,0)
 178       else if (index(controlcard,'CSA').gt.0) then
 179         modecalc=8
 180 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 181 crc
 182 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 183 crc
 184 crc        modecalc=9
 185 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 186 cfmc        modecalc=10
 187       else if (index(controlcard,'SOFTREG').gt.0) then
 188         modecalc=11
 189       else if (index(controlcard,'CHECK_BOND').gt.0) then
 190         modecalc=-1
 191       else if (index(controlcard,'TEST').gt.0) then
 192         modecalc=-2
 193       else if (index(controlcard,'MD').gt.0) then
 194         modecalc=12
 195       else if (index(controlcard,'RE ').gt.0) then
 196         modecalc=14
 197       endif
 198
 199       lmuca=index(controlcard,'MUCA').gt.0
 200       call readi(controlcard,'MUCADYN',mucadyn,0)
 201       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 202       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 203      & then
 204        write (iout,*) 'MUCADYN=',mucadyn
 205        write (iout,*) 'MUCASMOOTH=',muca_smooth
 206       endif
 207
 208       iscode=index(controlcard,'ONE_LETTER')
 209       indphi=index(controlcard,'PHI')
 210       indback=index(controlcard,'BACK')
 211       iranconf=index(controlcard,'RAND_CONF')
 212       i2ndstr=index(controlcard,'USE_SEC_PRED')
 213       gradout=index(controlcard,'GRADOUT').gt.0
 214       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 215
 216       if(me.eq.king.or..not.out1file)
 217      & write (iout,'(2a)') diagmeth(kdiag),
 218      &  ' routine used to diagonalize matrices.'
 219       return
 220       end
 221 c------------------------------------------------------------------------------
 222       subroutine molread
 223 C
 224 C Read molecular data.
 225 C
 226       implicit real*8 (a-h,o-z)
 227       include 'DIMENSIONS'
 228 #ifdef MPI
 229       include 'mpif.h'
 230       integer error_msg
 231 #endif
 232       include 'COMMON.IOUNITS'
 233       include 'COMMON.GEO'
 234       include 'COMMON.VAR'
 235       include 'COMMON.INTERACT'
 236       include 'COMMON.LOCAL'
 237       include 'COMMON.NAMES'
 238       include 'COMMON.CHAIN'
 239       include 'COMMON.FFIELD'
 240       include 'COMMON.SBRIDGE'
 241       include 'COMMON.HEADER'
 242       include 'COMMON.CONTROL'
 243 c      include 'COMMON.DBASE'
 244 c      include 'COMMON.THREAD'
 245       include 'COMMON.CONTACTS'
 246       include 'COMMON.TORCNSTR'
 247       include 'COMMON.TIME1'
 248       include 'COMMON.BOUNDS'
 249 c      include 'COMMON.MD'
 250 c      include 'COMMON.REMD'
 251       include 'COMMON.SETUP'
 252       character*4 sequence(maxres)
 253       integer rescode
 254       double precision x(maxvar)
 255       character*256 pdbfile
 256       character*320 weightcard
 257       character*80 weightcard_t,ucase
 258       dimension itype_pdb(maxres)
 259       common /pizda/ itype_pdb
 260       logical seq_comp,fail
 261       double precision energia(0:n_ene)
 262       integer ilen
 263       external ilen
 264 C
 265 C Body
 266 C
 267 C Read weights of the subsequent energy terms.
 268
 269        call card_concat(weightcard)
 270        call reada(weightcard,'WLONG',wlong,1.0D0)
 271        call reada(weightcard,'WSC',wsc,wlong)
 272        call reada(weightcard,'WSCP',wscp,wlong)
 273        call reada(weightcard,'WELEC',welec,1.0D0)
 274        call reada(weightcard,'WVDWPP',wvdwpp,welec)
 275        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 276        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 277        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 278        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 279        call reada(weightcard,'WTURN3',wturn3,1.0D0)
 280        call reada(weightcard,'WTURN4',wturn4,1.0D0)
 281        call reada(weightcard,'WTURN6',wturn6,1.0D0)
 282        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 283        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 284        call reada(weightcard,'WBOND',wbond,1.0D0)
 285        call reada(weightcard,'WTOR',wtor,1.0D0)
 286        call reada(weightcard,'WTORD',wtor_d,1.0D0)
 287        call reada(weightcard,'WANG',wang,1.0D0)
 288        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 289 C     Juyong
 290        call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
 291        call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
 292        call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
 293        call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
 294 C
 295        call reada(weightcard,'SCAL14',scal14,0.4D0)
 296        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 297        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 298        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 299        call reada(weightcard,'TEMP0',temp0,300.0d0)
 300        if (index(weightcard,'SOFT').gt.0) ipot=6
 301 C 12/1/95 Added weight for the multi-body term WCORR
 302        call reada(weightcard,'WCORRH',wcorr,1.0D0)
 303        if (wcorr4.gt.0.0d0) wcorr=wcorr4
 304        weights(1)=wsc
 305        weights(2)=wscp
 306        weights(3)=welec
 307        weights(4)=wcorr
 308        weights(5)=wcorr5
 309        weights(6)=wcorr6
 310        weights(7)=wel_loc
 311        weights(8)=wturn3
 312        weights(9)=wturn4
 313        weights(10)=wturn6
 314        weights(11)=wang
 315        weights(12)=wscloc
 316        weights(13)=wtor
 317        weights(14)=wtor_d
 318        weights(15)=wstrain
 319        weights(16)=wvdwpp
 320        weights(17)=wbond
 321        weights(18)=scal14
 322        weights(21)=wsccor
 323 C     JUYONG
 324        weights(24)=wdfa_dist
 325        weights(25)=wdfa_tor
 326        weights(26)=wdfa_nei
 327        weights(27)=wdfa_beta
 328 C
 329
 330       if(me.eq.king.or..not.out1file)
 331      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 332      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 333      &  wturn4,wturn6,
 334      &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
 335
 336    10 format (/'Energy-term weights (unscaled):'//
 337      & 'WSCC=   ',f10.6,' (SC-SC)'/
 338      & 'WSCP=   ',f10.6,' (SC-p)'/
 339      & 'WELEC=  ',f10.6,' (p-p electr)'/
 340      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 341      & 'WBOND=  ',f10.6,' (stretching)'/
 342      & 'WANG=   ',f10.6,' (bending)'/
 343      & 'WSCLOC= ',f10.6,' (SC local)'/
 344      & 'WTOR=   ',f10.6,' (torsional)'/
 345      & 'WTORD=  ',f10.6,' (double torsional)'/
 346      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 347      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 348      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 349      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 350      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 351      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
 352      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 353      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 354      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 355      & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
 356      & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
 357      & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
 358      & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
 359
 360       if(me.eq.king.or..not.out1file)then
 361        if (wcorr4.gt.0.0d0) then
 362         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 363      &   'between contact pairs of peptide groups'
 364         write (iout,'(2(a,f5.3/))')
 365      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 366      &  'Range of quenching the correlation terms:',2*delt_corr
 367        else if (wcorr.gt.0.0d0) then
 368         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 369      &   'between contact pairs of peptide groups'
 370        endif
 371        write (iout,'(a,f8.3)')
 372      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 373        write (iout,'(a,f8.3)')
 374      &  'General scaling factor of SC-p interactions:',scalscp
 375       endif
 376       r0_corr=cutoff_corr-delt_corr
 377       do i=1,20
 378         aad(i,1)=scalscp*aad(i,1)
 379         aad(i,2)=scalscp*aad(i,2)
 380         bad(i,1)=scalscp*bad(i,1)
 381         bad(i,2)=scalscp*bad(i,2)
 382       enddo
 383 c      call rescale_weights(t_bath)
 384       if(me.eq.king.or..not.out1file)
 385      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 386      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 387      &  wturn4,wturn6,
 388      &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
 389
 390    22 format (/'Energy-term weights (scaled):'//
 391      & 'WSCC=   ',f10.6,' (SC-SC)'/
 392      & 'WSCP=   ',f10.6,' (SC-p)'/
 393      & 'WELEC=  ',f10.6,' (p-p electr)'/
 394      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 395      & 'WBOND=  ',f10.6,' (stretching)'/
 396      & 'WANG=   ',f10.6,' (bending)'/
 397      & 'WSCLOC= ',f10.6,' (SC local)'/
 398      & 'WTOR=   ',f10.6,' (torsional)'/
 399      & 'WTORD=  ',f10.6,' (double torsional)'/
 400      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 401      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 402      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 403      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 404      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 405      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
 406      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 407      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 408      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 409      & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
 410      & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
 411      & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
 412      & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
 413
 414       if(me.eq.king.or..not.out1file)
 415      & write (iout,*) "Reference temperature for weights calculation:",
 416      &  temp0
 417       call reada(weightcard,"D0CM",d0cm,3.78d0)
 418       call reada(weightcard,"AKCM",akcm,15.1d0)
 419       call reada(weightcard,"AKTH",akth,11.0d0)
 420       call reada(weightcard,"AKCT",akct,12.0d0)
 421       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 422       call reada(weightcard,"V2SS",v2ss,7.61d0)
 423       call reada(weightcard,"V3SS",v3ss,13.7d0)
 424       call reada(weightcard,"EBR",ebr,-5.50D0)
 425       if(me.eq.king.or..not.out1file) then
 426        write (iout,*) "Parameters of the SS-bond potential:"
 427        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 428      & " AKCT",akct
 429        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 430        write (iout,*) "EBR",ebr
 431        print *,'indpdb=',indpdb,' pdbref=',pdbref
 432       endif
 433       if (indpdb.gt.0 .or. pdbref) then
 434         read(inp,'(a)') pdbfile
 435         if(me.eq.king.or..not.out1file)
 436      &   write (iout,'(2a)') 'PDB data will be read from file ',
 437      &   pdbfile(:ilen(pdbfile))
 438         open(ipdbin,file=pdbfile,status='old',err=33)
 439         goto 34
 440   33    write (iout,'(a)') 'Error opening PDB file.'
 441         stop
 442   34    continue
 443 c        print *,'Begin reading pdb data'
 444         call readpdb
 445 c        print *,'Finished reading pdb data'
 446         if(me.eq.king.or..not.out1file)
 447      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 448      &   ' nstart_sup=',nstart_sup
 449         do i=1,nres
 450           itype_pdb(i)=itype(i)
 451         enddo
 452         close (ipdbin)
 453         nnt=nstart_sup
 454         nct=nstart_sup+nsup-1
 455         call contact(.false.,ncont_ref,icont_ref,co)
 456
 457         if (sideadd) then
 458 C Following 2 lines for diagnostics; comment out if not needed
 459          write (iout,*) "Before sideadd"
 460          call intout
 461          if(me.eq.king.or..not.out1file)
 462      &    write(iout,*)'Adding sidechains'
 463          maxsi=1000
 464          do i=2,nres-1
 465           iti=itype(i)
 466           if (iti.ne.10) then
 467             nsi=0
 468             fail=.true.
 469             do while (fail.and.nsi.le.maxsi)
 470 c              call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 471               nsi=nsi+1
 472             enddo
 473             if(fail) write(iout,*)'Adding sidechain failed for res ',
 474      &              i,' after ',nsi,' trials'
 475           endif
 476          enddo
 477         endif
 478 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
 479         call chainbuild
 480 C Following 2 lines for diagnostics; comment out if not needed
 481         write (iout,*) "After sideadd"
 482         call intout
 483       endif
 484
 485       if (indpdb.eq.0) then
 486 C Read sequence if not taken from the pdb file.
 487         read (inp,*) nres
 488 c        print *,'nres=',nres
 489         if (iscode.gt.0) then
 490           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 491         else
 492           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 493         endif
 494 C Convert sequence to numeric code
 495         do i=1,nres
 496           itype(i)=rescode(i,sequence(i),iscode)
 497         enddo
 498 C Assign initial virtual bond lengths
 499         do i=2,nres
 500           vbld(i)=vbl
 501           vbld_inv(i)=vblinv
 502         enddo
 503         do i=2,nres-1
 504           vbld(i+nres)=dsc(itype(i))
 505           vbld_inv(i+nres)=dsc_inv(itype(i))
 506 c          write (iout,*) "i",i," itype",itype(i),
 507 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 508         enddo
 509       endif
 510 c      print *,nres
 511 c      print '(20i4)',(itype(i),i=1,nres)
 512       do i=1,nres
 513 #ifdef PROCOR
 514         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 515 #else
 516         if (itype(i).eq.21) then
 517 #endif
 518           itel(i)=0
 519 #ifdef PROCOR
 520         else if (itype(i+1).ne.20) then
 521 #else
 522         else if (itype(i).ne.20) then
 523 #endif
 524           itel(i)=1
 525         else
 526           itel(i)=2
 527         endif
 528       enddo
 529       if(me.eq.king.or..not.out1file)then
 530        write (iout,*) "ITEL"
 531        do i=1,nres-1
 532          write (iout,*) i,itype(i),itel(i)
 533        enddo
 534        print *,'Call Read_Bridge.'
 535       endif
 536       call read_bridge
 537 C 8/13/98 Set limits to generating the dihedral angles
 538       do i=1,nres
 539         phibound(1,i)=-pi
 540         phibound(2,i)=pi
 541       enddo
 542       read (inp,*) ndih_constr
 543       if (ndih_constr.gt.0) then
 544         read (inp,*) ftors
 545         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 546         if(me.eq.king.or..not.out1file)then
 547          write (iout,*)
 548      &   'There are',ndih_constr,' constraints on phi angles.'
 549          do i=1,ndih_constr
 550           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 551          enddo
 552         endif
 553         do i=1,ndih_constr
 554           phi0(i)=deg2rad*phi0(i)
 555           drange(i)=deg2rad*drange(i)
 556         enddo
 557         if(me.eq.king.or..not.out1file)
 558      &   write (iout,*) 'FTORS',ftors
 559         do i=1,ndih_constr
 560           ii = idih_constr(i)
 561           phibound(1,ii) = phi0(i)-drange(i)
 562           phibound(2,ii) = phi0(i)+drange(i)
 563         enddo
 564       endif
 565       nnt=1
 566 #ifdef MPI
 567       if (me.eq.king) then
 568 #endif
 569        write (iout,'(a)') 'Boundaries in phi angle sampling:'
 570        do i=1,nres
 571          write (iout,'(a3,i5,2f10.1)')
 572      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 573        enddo
 574 #ifdef MP
 575       endif
 576 #endif
 577       nct=nres
 578 cd      print *,'NNT=',NNT,' NCT=',NCT
 579       if (itype(1).eq.21) nnt=2
 580       if (itype(nres).eq.21) nct=nct-1
 581
 582 C     Juyong:READ init_vars
 583 C     Initialize variables!
 584 C     Juyong:READ read_info
 585 C     READ fragment information!!
 586 C     both routines should be in dfa.F file!!
 587
 588       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
 589      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
 590        call init_dfa_vars
 591        print*, 'init_dfa_vars finished!'
 592        call read_dfa_info
 593        print*, 'read_dfa_info finished!'
 594       endif
 595 C
 596 C
 597
 598
 599       if (pdbref) then
 600         if(me.eq.king.or..not.out1file)
 601      &   write (iout,'(a,i3)') 'nsup=',nsup
 602         nstart_seq=nnt
 603         if (nsup.le.(nct-nnt+1)) then
 604           do i=0,nct-nnt+1-nsup
 605             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 606               nstart_seq=nnt+i
 607               goto 111
 608             endif
 609           enddo
 610           write (iout,'(a)')
 611      &            'Error - sequences to be superposed do not match.'
 612           stop
 613         else
 614           do i=0,nsup-(nct-nnt+1)
 615             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 616      &      then
 617               nstart_sup=nstart_sup+i
 618               nsup=nct-nnt+1
 619               goto 111
 620             endif
 621           enddo
 622           write (iout,'(a)')
 623      &            'Error - sequences to be superposed do not match.'
 624         endif
 625   111   continue
 626         if (nsup.eq.0) nsup=nct-nnt
 627         if (nstart_sup.eq.0) nstart_sup=nnt
 628         if (nstart_seq.eq.0) nstart_seq=nnt
 629         if(me.eq.king.or..not.out1file)
 630      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 631      &                 ' nstart_seq=',nstart_seq
 632       endif
 633 c--- Zscore rms -------
 634       if (nz_start.eq.0) nz_start=nnt
 635       if (nz_end.eq.0 .and. nsup.gt.0) then
 636         nz_end=nnt+nsup-1
 637       else if (nz_end.eq.0) then
 638         nz_end=nct
 639       endif
 640       if(me.eq.king.or..not.out1file)then
 641        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
 642        write (iout,*) 'IZ_SC=',iz_sc
 643       endif
 644 c----------------------
 645       call init_int_table
 646       if (refstr) then
 647         if (.not.pdbref) then
 648           call read_angles(inp,*38)
 649           goto 39
 650    38     write (iout,'(a)') 'Error reading reference structure.'
 651 #ifdef MPI
 652           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 653           stop 'Error reading reference structure'
 654 #endif
 655    39     call chainbuild
 656           call setup_var
 657 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
 658           nstart_sup=nnt
 659           nstart_seq=nnt
 660           nsup=nct-nnt+1
 661           do i=1,2*nres
 662             do j=1,3
 663               cref(j,i)=c(j,i)
 664             enddo
 665           enddo
 666           call contact(.true.,ncont_ref,icont_ref,co)
 667         endif
 668         if(me.eq.king.or..not.out1file)
 669      &   write (iout,*) 'Contact order:',co
 670         if (pdbref) then
 671         if(me.eq.king.or..not.out1file)
 672      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
 673         do i=1,ncont_ref
 674           do j=1,2
 675             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
 676           enddo
 677           if(me.eq.king.or..not.out1file)
 678      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
 679      &     icont_ref(1,i),' ',
 680      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
 681         enddo
 682         endif
 683       endif
 684 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
 685       if (constr_dist.gt.0) then
 686         call read_dist_constr
 687       endif
 688       if (nhpb.gt.0) call hpb_partition
 689 c      write (iout,*) "After read_dist_constr nhpb",nhpb
 690 c      call flush(iout)
 691       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
 692      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
 693      &    modecalc.ne.10) then
 694 C If input structure hasn't been supplied from the PDB file read or generate
 695 C initial geometry.
 696         if (iranconf.eq.0 .and. .not. extconf) then
 697           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
 698      &     write (iout,'(a)') 'Initial geometry will be read in.'
 699           if (read_cart) then
 700             read(inp,'(8f10.5)',end=36,err=36)
 701      &       ((c(l,k),l=1,3),k=1,nres),
 702      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 703             call int_from_cart1(.false.)
 704             do i=1,nres-1
 705               do j=1,3
 706                 dc(j,i)=c(j,i+1)-c(j,i)
 707                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
 708               enddo
 709             enddo
 710             do i=nnt,nct
 711               if (itype(i).ne.10) then
 712                 do j=1,3
 713                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
 714                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
 715                 enddo
 716               endif
 717             enddo
 718             return
 719           else
 720             call read_angles(inp,*36)
 721           endif
 722           goto 37
 723    36     write (iout,'(a)') 'Error reading angle file.'
 724 #ifdef MPI
 725           call mpi_finalize( MPI_COMM_WORLD,IERR )
 726 #endif
 727           stop 'Error reading angle file.'
 728    37     continue
 729         else if (extconf) then
 730          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
 731      &    write (iout,'(a)') 'Extended chain initial geometry.'
 732          do i=3,nres
 733           theta(i)=90d0*deg2rad
 734          enddo
 735          do i=4,nres
 736           phi(i)=180d0*deg2rad
 737          enddo
 738          do i=2,nres-1
 739           alph(i)=110d0*deg2rad
 740          enddo
 741          do i=2,nres-1
 742           omeg(i)=-120d0*deg2rad
 743          enddo
 744         else
 745           if(me.eq.king.or..not.out1file)
 746      &     write (iout,'(a)') 'Random-generated initial geometry.'
 747
 748
 749 #ifdef MPI
 750           if (me.eq.king  .or. fg_rank.eq.0 .and. (
 751      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
 752 #endif
 753             do itrial=1,100
 754               itmp=1
 755 c              call gen_rand_conf(itmp,*30)
 756               goto 40
 757    30         write (iout,*) 'Failed to generate random conformation',
 758      &          ', itrial=',itrial
 759               write (*,*) 'Processor:',me,
 760      &          ' Failed to generate random conformation',
 761      &          ' itrial=',itrial
 762               call intout
 763
 764 #ifdef AIX
 765               call flush_(iout)
 766 #else
 767               call flush(iout)
 768 #endif
 769             enddo
 770             write (iout,'(a,i3,a)') 'Processor:',me,
 771      &        ' error in generating random conformation.'
 772             write (*,'(a,i3,a)') 'Processor:',me,
 773      &        ' error in generating random conformation.'
 774             call flush(iout)
 775 #ifdef MPI
 776             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
 777    40       continue
 778           endif
 779 #else
 780           do itrial=1,100
 781             itmp=1
 782 c            call gen_rand_conf(itmp,*31)
 783             goto 40
 784    31       write (iout,*) 'Failed to generate random conformation',
 785      &        ', itrial=',itrial
 786             write (*,*) 'Failed to generate random conformation',
 787      &        ', itrial=',itrial
 788           enddo
 789           write (iout,'(a,i3,a)') 'Processor:',me,
 790      &      ' error in generating random conformation.'
 791           write (*,'(a,i3,a)') 'Processor:',me,
 792      &      ' error in generating random conformation.'
 793           stop
 794    40     continue
 795 #endif
 796         endif
 797       elseif (modecalc.eq.4) then
 798         read (inp,'(a)') intinname
 799         open (intin,file=intinname,status='old',err=333)
 800         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
 801      &  write (iout,'(a)') 'intinname',intinname
 802         write (*,'(a)') 'Processor',myrank,' intinname',intinname
 803         goto 334
 804   333   write (iout,'(2a)') 'Error opening angle file ',intinname
 805 #ifdef MPI
 806         call MPI_Finalize(MPI_COMM_WORLD,IERR)
 807 #endif
 808         stop 'Error opening angle file.'
 809   334   continue
 810
 811       endif
 812 C Generate distance constraints, if the PDB structure is to be regularized.
 813       call setup_var
 814       if (me.eq.king .or. .not. out1file)
 815      & call intout
 816       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
 817         write (iout,'(/a,i3,a)')
 818      &  'The chain contains',ns,' disulfide-bridging cysteines.'
 819         write (iout,'(20i4)') (iss(i),i=1,ns)
 820         write (iout,'(/a/)') 'Pre-formed links are:'
 821         do i=1,nss
 822           i1=ihpb(i)-nres
 823           i2=jhpb(i)-nres
 824           it1=itype(i1)
 825           it2=itype(i2)
 826           if (me.eq.king.or..not.out1file)
 827      &    write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 828      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
 829      &    ebr,forcon(i)
 830         enddo
 831         write (iout,'(a)')
 832       endif
 833 c       if (i2ndstr.gt.0) call secstrp2dihc
 834 c      call geom_to_var(nvar,x)
 835 c      call etotal(energia(0))
 836 c      call enerprint(energia(0))
 837 c      call briefout(0,etot)
 838 c      stop
 839 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
 840 cd    write (iout,'(a)') 'Variable list:'
 841 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
 842 #ifdef MPI
 843       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
 844      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
 845      &  'Processor',myrank,': end reading molecular data.'
 846 #endif
 847       return
 848       end
 849 c--------------------------------------------------------------------------
 850       logical function seq_comp(itypea,itypeb,length)
 851       implicit none
 852       integer length,itypea(length),itypeb(length)
 853       integer i
 854       do i=1,length
 855         if (itypea(i).ne.itypeb(i)) then
 856           seq_comp=.false.
 857           return
 858         endif
 859       enddo
 860       seq_comp=.true.
 861       return
 862       end
 863 c-----------------------------------------------------------------------------
 864       subroutine read_bridge
 865 C Read information about disulfide bridges.
 866       implicit real*8 (a-h,o-z)
 867       include 'DIMENSIONS'
 868 #ifdef MPI
 869       include 'mpif.h'
 870 #endif
 871       include 'COMMON.IOUNITS'
 872       include 'COMMON.GEO'
 873       include 'COMMON.VAR'
 874       include 'COMMON.INTERACT'
 875       include 'COMMON.LOCAL'
 876       include 'COMMON.NAMES'
 877       include 'COMMON.CHAIN'
 878       include 'COMMON.FFIELD'
 879       include 'COMMON.SBRIDGE'
 880       include 'COMMON.HEADER'
 881       include 'COMMON.CONTROL'
 882 c      include 'COMMON.DBASE'
 883 c      include 'COMMON.THREAD'
 884       include 'COMMON.TIME1'
 885       include 'COMMON.SETUP'
 886 C Read bridging residues.
 887       read (inp,*) ns,(iss(i),i=1,ns)
 888       print *,'ns=',ns
 889       if(me.eq.king.or..not.out1file)
 890      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 891 C Check whether the specified bridging residues are cystines.
 892       do i=1,ns
 893         if (itype(iss(i)).ne.1) then
 894           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
 895      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 896      &   ' can form a disulfide bridge?!!!'
 897           write (*,'(2a,i3,a)')
 898      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 899      &   ' can form a disulfide bridge?!!!'
 900 #ifdef MPI
 901          call MPI_Finalize(MPI_COMM_WORLD,ierror)
 902          stop
 903 #endif
 904         endif
 905       enddo
 906 C Read preformed bridges.
 907       if (ns.gt.0) then
 908       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 909       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 910       if (nss.gt.0) then
 911         nhpb=nss
 912 C Check if the residues involved in bridges are in the specified list of
 913 C bridging residues.
 914         do i=1,nss
 915           do j=1,i-1
 916             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 917      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 918               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 919      &      ' contains residues present in other pairs.'
 920               write (*,'(a,i3,a)') 'Disulfide pair',i,
 921      &      ' contains residues present in other pairs.'
 922 #ifdef MPI
 923               call MPI_Finalize(MPI_COMM_WORLD,ierror)
 924               stop
 925 #endif
 926             endif
 927           enddo
 928           do j=1,ns
 929             if (ihpb(i).eq.iss(j)) goto 10
 930           enddo
 931           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 932    10     continue
 933           do j=1,ns
 934             if (jhpb(i).eq.iss(j)) goto 20
 935           enddo
 936           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 937    20     continue
 938           dhpb(i)=dbr
 939           forcon(i)=fbr
 940         enddo
 941         do i=1,nss
 942           ihpb(i)=ihpb(i)+nres
 943           jhpb(i)=jhpb(i)+nres
 944         enddo
 945       endif
 946       endif
 947       return
 948       end
 949 c----------------------------------------------------------------------------
 950       subroutine read_x(kanal,*)
 951       implicit real*8 (a-h,o-z)
 952       include 'DIMENSIONS'
 953       include 'COMMON.GEO'
 954       include 'COMMON.VAR'
 955       include 'COMMON.CHAIN'
 956       include 'COMMON.IOUNITS'
 957       include 'COMMON.CONTROL'
 958       include 'COMMON.LOCAL'
 959       include 'COMMON.INTERACT'
 960 c Read coordinates from input
 961 c
 962       read(kanal,'(8f10.5)',end=10,err=10)
 963      &  ((c(l,k),l=1,3),k=1,nres),
 964      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
 965       do j=1,3
 966         c(j,nres+1)=c(j,1)
 967         c(j,2*nres)=c(j,nres)
 968       enddo
 969       call int_from_cart1(.false.)
 970       do i=1,nres-1
 971         do j=1,3
 972           dc(j,i)=c(j,i+1)-c(j,i)
 973           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 974         enddo
 975       enddo
 976       do i=nnt,nct
 977         if (itype(i).ne.10) then
 978           do j=1,3
 979             dc(j,i+nres)=c(j,i+nres)-c(j,i)
 980             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 981           enddo
 982         endif
 983       enddo
 984
 985       return
 986    10 return1
 987       end
 988 c------------------------------------------------------------------------------
 989       subroutine setup_var
 990       implicit real*8 (a-h,o-z)
 991       include 'DIMENSIONS'
 992       include 'COMMON.IOUNITS'
 993       include 'COMMON.GEO'
 994       include 'COMMON.VAR'
 995       include 'COMMON.INTERACT'
 996       include 'COMMON.LOCAL'
 997       include 'COMMON.NAMES'
 998       include 'COMMON.CHAIN'
 999       include 'COMMON.FFIELD'
1000       include 'COMMON.SBRIDGE'
1001       include 'COMMON.HEADER'
1002       include 'COMMON.CONTROL'
1003 c      include 'COMMON.DBASE'
1004 c      include 'COMMON.THREAD'
1005       include 'COMMON.TIME1'
1006 C Set up variable list.
1007       ntheta=nres-2
1008       nphi=nres-3
1009       nvar=ntheta+nphi
1010       nside=0
1011       do i=2,nres-1
1012         if (itype(i).ne.10) then
1013           nside=nside+1
1014           ialph(i,1)=nvar+nside
1015           ialph(nside,2)=i
1016         endif
1017       enddo
1018       if (indphi.gt.0) then
1019         nvar=nphi
1020       else if (indback.gt.0) then
1021         nvar=nphi+ntheta
1022       else
1023         nvar=nvar+2*nside
1024       endif
1025 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1026       return
1027       end
1028 c----------------------------------------------------------------------------
1029       subroutine gen_dist_constr
1030 C Generate CA distance constraints.
1031       implicit real*8 (a-h,o-z)
1032       include 'DIMENSIONS'
1033       include 'COMMON.IOUNITS'
1034       include 'COMMON.GEO'
1035       include 'COMMON.VAR'
1036       include 'COMMON.INTERACT'
1037       include 'COMMON.LOCAL'
1038       include 'COMMON.NAMES'
1039       include 'COMMON.CHAIN'
1040       include 'COMMON.FFIELD'
1041       include 'COMMON.SBRIDGE'
1042       include 'COMMON.HEADER'
1043       include 'COMMON.CONTROL'
1044 c      include 'COMMON.DBASE'
1045 c      include 'COMMON.THREAD'
1046       include 'COMMON.TIME1'
1047       dimension itype_pdb(maxres)
1048       common /pizda/ itype_pdb
1049       character*2 iden
1050 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1051 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1052 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1053 cd     & ' nsup',nsup
1054       do i=nstart_sup,nstart_sup+nsup-1
1055 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1056 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1057         do j=i+2,nstart_sup+nsup-1
1058           nhpb=nhpb+1
1059           ihpb(nhpb)=i+nstart_seq-nstart_sup
1060           jhpb(nhpb)=j+nstart_seq-nstart_sup
1061           forcon(nhpb)=weidis
1062           dhpb(nhpb)=dist(i,j)
1063         enddo
1064       enddo
1065 cd      write (iout,'(a)') 'Distance constraints:'
1066 cd      do i=nss+1,nhpb
1067 cd        ii=ihpb(i)
1068 cd        jj=jhpb(i)
1069 cd        iden='CA'
1070 cd        if (ii.gt.nres) then
1071 cd          iden='SC'
1072 cd          ii=ii-nres
1073 cd          jj=jj-nres
1074 cd        endif
1075 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1076 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1077 cd     &  dhpb(i),forcon(i)
1078 cd      enddo
1079       return
1080       end
1081 c----------------------------------------------------------------------------
1082       subroutine csaread
1083       implicit real*8 (a-h,o-z)
1084       include 'DIMENSIONS'
1085       include 'COMMON.IOUNITS'
1086       include 'COMMON.GEO'
1087       include 'COMMON.CSA'
1088       include 'COMMON.BANK'
1089       include 'COMMON.CONTROL'
1090       character*80 ucase
1091       character*620 mcmcard
1092       call card_concat(mcmcard)
1093
1094       call readi(mcmcard,'NCONF',nconf,50)
1095       call readi(mcmcard,'NADD',nadd,0)
1096       call readi(mcmcard,'JSTART',jstart,1)
1097       call readi(mcmcard,'JEND',jend,1)
1098       call readi(mcmcard,'NSTMAX',nstmax,500000)
1099       call readi(mcmcard,'N0',n0,1)
1100       call readi(mcmcard,'N1',n1,6)
1101       call readi(mcmcard,'N2',n2,4)
1102       call readi(mcmcard,'N3',n3,0)
1103       call readi(mcmcard,'N4',n4,0)
1104       call readi(mcmcard,'N5',n5,0)
1105       call readi(mcmcard,'N6',n6,10)
1106       call readi(mcmcard,'N7',n7,0)
1107       call readi(mcmcard,'N8',n8,0)
1108       call readi(mcmcard,'N9',n9,0)
1109       call readi(mcmcard,'N14',n14,0)
1110       call readi(mcmcard,'N15',n15,0)
1111       call readi(mcmcard,'N16',n16,0)
1112       call readi(mcmcard,'N17',n17,0)
1113       call readi(mcmcard,'N18',n18,0)
1114
1115       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1116
1117       call readi(mcmcard,'NDIFF',ndiff,2)
1118       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1119       call readi(mcmcard,'IS1',is1,1)
1120       call readi(mcmcard,'IS2',is2,8)
1121       call readi(mcmcard,'NRAN0',nran0,4)
1122       call readi(mcmcard,'NRAN1',nran1,2)
1123       call readi(mcmcard,'IRR',irr,1)
1124       call readi(mcmcard,'NSEED',nseed,20)
1125       call readi(mcmcard,'NTOTAL',ntotal,10000)
1126       call reada(mcmcard,'CUT1',cut1,2.0d0)
1127       call reada(mcmcard,'CUT2',cut2,5.0d0)
1128       call reada(mcmcard,'ESTOP',estop,-300000.0d0)
1129       call readi(mcmcard,'ICMAX',icmax,1)
1130       call readi(mcmcard,'NBANKM',nbankm,400)
1131       call readi(mcmcard,'IUCUT',iucut,2)
1132       call readi(mcmcard,'IRESTART',irestart,0)
1133 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1134       ntbankm=0
1135 c!bankt
1136       call reada(mcmcard,'DELE',dele,20.0d0)
1137       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1138       call readi(mcmcard,'IREF',iref,0)
1139       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1140       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1141       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1142       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1143       write (iout,*) "NCONF_IN",nconf_in
1144       tm_score=(index(mcmcard,'TMSCORE').gt.0)
1145       if (tm_score) write (iout,*) "Using TM_Score instead of DIFF",
1146      &      " for torsional angles"
1147       return
1148       end
1149
1150       subroutine read_minim
1151       implicit real*8 (a-h,o-z)
1152       include 'DIMENSIONS'
1153       include 'COMMON.MINIM'
1154       include 'COMMON.IOUNITS'
1155       character*80 ucase
1156       character*320 minimcard
1157       call card_concat(minimcard)
1158       call readi(minimcard,'MAXMIN',maxmin,2000)
1159       call readi(minimcard,'MAXFUN',maxfun,5000)
1160       call readi(minimcard,'MINMIN',minmin,maxmin)
1161       call readi(minimcard,'MINFUN',minfun,maxmin)
1162       call reada(minimcard,'TOLF',tolf,1.0D-2)
1163       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1164       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1165       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1166       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1167       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1168      &         'Options in energy minimization:'
1169       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1170      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1171      & 'MinMin:',MinMin,' MinFun:',MinFun,
1172      & ' TolF:',TolF,' RTolF:',RTolF
1173       return
1174       end
1175 c----------------------------------------------------------------------------
1176       subroutine read_angles(kanal,*)
1177       implicit real*8 (a-h,o-z)
1178       include 'DIMENSIONS'
1179       include 'COMMON.GEO'
1180       include 'COMMON.VAR'
1181       include 'COMMON.CHAIN'
1182       include 'COMMON.IOUNITS'
1183       include 'COMMON.CONTROL'
1184 c Read angles from input
1185 c
1186        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1187        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1188        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1189        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1190
1191        do i=1,nres
1192 c 9/7/01 avoid 180 deg valence angle
1193         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1194 c
1195         theta(i)=deg2rad*theta(i)
1196         phi(i)=deg2rad*phi(i)
1197         alph(i)=deg2rad*alph(i)
1198         omeg(i)=deg2rad*omeg(i)
1199        enddo
1200       return
1201    10 return1
1202       end
1203 c----------------------------------------------------------------------------
1204       subroutine reada(rekord,lancuch,wartosc,default)
1205       implicit none
1206       character*(*) rekord,lancuch
1207       double precision wartosc,default
1208       integer ilen,iread
1209       external ilen
1210       iread=index(rekord,lancuch)
1211       if (iread.eq.0) then
1212         wartosc=default
1213         return
1214       endif
1215       iread=iread+ilen(lancuch)+1
1216       read (rekord(iread:),*,err=10,end=10) wartosc
1217       return
1218   10  wartosc=default
1219       return
1220       end
1221 c----------------------------------------------------------------------------
1222       subroutine readi(rekord,lancuch,wartosc,default)
1223       implicit none
1224       character*(*) rekord,lancuch
1225       integer wartosc,default
1226       integer ilen,iread
1227       external ilen
1228       iread=index(rekord,lancuch)
1229       if (iread.eq.0) then
1230         wartosc=default
1231         return
1232       endif
1233       iread=iread+ilen(lancuch)+1
1234       read (rekord(iread:),*,err=10,end=10) wartosc
1235       return
1236   10  wartosc=default
1237       return
1238       end
1239 c----------------------------------------------------------------------------
1240       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1241       implicit none
1242       integer dim,i
1243       integer tablica(dim),default
1244       character*(*) rekord,lancuch
1245       character*80 aux
1246       integer ilen,iread
1247       external ilen
1248       do i=1,dim
1249         tablica(i)=default
1250       enddo
1251       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1252       if (iread.eq.0) return
1253       iread=iread+ilen(lancuch)+1
1254       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1255    10 return
1256       end
1257 c----------------------------------------------------------------------------
1258       subroutine multreada(rekord,lancuch,tablica,dim,default)
1259       implicit none
1260       integer dim,i
1261       double precision tablica(dim),default
1262       character*(*) rekord,lancuch
1263       character*80 aux
1264       integer ilen,iread
1265       external ilen
1266       do i=1,dim
1267         tablica(i)=default
1268       enddo
1269       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1270       if (iread.eq.0) return
1271       iread=iread+ilen(lancuch)+1
1272       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1273    10 return
1274       end
1275 c----------------------------------------------------------------------------
1276       subroutine openunits
1277       implicit real*8 (a-h,o-z)
1278       include 'DIMENSIONS'
1279 #ifdef MPI
1280       include 'mpif.h'
1281       character*16 form,nodename
1282       integer nodelen
1283 #endif
1284       include 'COMMON.SETUP'
1285       include 'COMMON.IOUNITS'
1286 c      include 'COMMON.MD'
1287       include 'COMMON.CONTROL'
1288       integer lenpre,lenpot,ilen,lentmp
1289       external ilen
1290       character*3 out1file_text,ucase
1291       character*3 ll
1292       external ucase
1293 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1294       call getenv_loc("PREFIX",prefix)
1295       pref_orig = prefix
1296       call getenv_loc("POT",pot)
1297       call getenv_loc("DIRTMP",tmpdir)
1298       call getenv_loc("CURDIR",curdir)
1299       call getenv_loc("OUT1FILE",out1file_text)
1300 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1301       out1file_text=ucase(out1file_text)
1302       if (out1file_text(1:1).eq."Y") then
1303         out1file=.true.
1304       else
1305         out1file=fg_rank.gt.0
1306       endif
1307       lenpre=ilen(prefix)
1308       lenpot=ilen(pot)
1309       lentmp=ilen(tmpdir)
1310       if (lentmp.gt.0) then
1311           write (*,'(80(1h!))')
1312           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1313           write (*,'(80(1h!))')
1314           write (*,*)"All output files will be on node /tmp directory."
1315 #ifdef MPI
1316         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1317         if (me.eq.king) then
1318           write (*,*) "The master node is ",nodename
1319         else if (fg_rank.eq.0) then
1320           write (*,*) "I am the CG slave node ",nodename
1321         else
1322           write (*,*) "I am the FG slave node ",nodename
1323         endif
1324 #endif
1325         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1326         lenpre = lentmp+lenpre+1
1327       endif
1328       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1329 C Get the names and open the input files
1330 #if defined(WINIFL) || defined(WINPGI)
1331       open(1,file=pref_orig(:ilen(pref_orig))//
1332      &  '.inp',status='old',readonly,shared)
1333        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1334 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1335 C Get parameter filenames and open the parameter files.
1336       call getenv_loc('BONDPAR',bondname)
1337       open (ibond,file=bondname,status='old',readonly,shared)
1338       call getenv_loc('THETPAR',thetname)
1339       open (ithep,file=thetname,status='old',readonly,shared)
1340 #ifndef CRYST_THETA
1341       call getenv_loc('THETPARPDB',thetname_pdb)
1342       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1343 #endif
1344       call getenv_loc('ROTPAR',rotname)
1345       open (irotam,file=rotname,status='old',readonly,shared)
1346 #ifndef CRYST_SC
1347       call getenv_loc('ROTPARPDB',rotname_pdb)
1348       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1349 #endif
1350       call getenv_loc('TORPAR',torname)
1351       open (itorp,file=torname,status='old',readonly,shared)
1352       call getenv_loc('TORDPAR',tordname)
1353       open (itordp,file=tordname,status='old',readonly,shared)
1354       call getenv_loc('FOURIER',fouriername)
1355       open (ifourier,file=fouriername,status='old',readonly,shared)
1356       call getenv_loc('ELEPAR',elename)
1357       open (ielep,file=elename,status='old',readonly,shared)
1358       call getenv_loc('SIDEPAR',sidename)
1359       open (isidep,file=sidename,status='old',readonly,shared)
1360 #elif (defined CRAY) || (defined AIX)
1361       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1362      &  action='read')
1363 c      print *,"Processor",myrank," opened file 1"
1364       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1365 c      print *,"Processor",myrank," opened file 9"
1366 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1367 C Get parameter filenames and open the parameter files.
1368       call getenv_loc('BONDPAR',bondname)
1369       open (ibond,file=bondname,status='old',action='read')
1370 c      print *,"Processor",myrank," opened file IBOND"
1371       call getenv_loc('THETPAR',thetname)
1372       open (ithep,file=thetname,status='old',action='read')
1373 c      print *,"Processor",myrank," opened file ITHEP"
1374 #ifndef CRYST_THETA
1375       call getenv_loc('THETPARPDB',thetname_pdb)
1376       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1377 #endif
1378       call getenv_loc('ROTPAR',rotname)
1379       open (irotam,file=rotname,status='old',action='read')
1380 c      print *,"Processor",myrank," opened file IROTAM"
1381 #ifndef CRYST_SC
1382       call getenv_loc('ROTPARPDB',rotname_pdb)
1383       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1384 #endif
1385       call getenv_loc('TORPAR',torname)
1386       open (itorp,file=torname,status='old',action='read')
1387 c      print *,"Processor",myrank," opened file ITORP"
1388       call getenv_loc('TORDPAR',tordname)
1389       open (itordp,file=tordname,status='old',action='read')
1390 c      print *,"Processor",myrank," opened file ITORDP"
1391       call getenv_loc('SCCORPAR',sccorname)
1392       open (isccor,file=sccorname,status='old',action='read')
1393 c      print *,"Processor",myrank," opened file ISCCOR"
1394       call getenv_loc('FOURIER',fouriername)
1395       open (ifourier,file=fouriername,status='old',action='read')
1396 c      print *,"Processor",myrank," opened file IFOURIER"
1397       call getenv_loc('ELEPAR',elename)
1398       open (ielep,file=elename,status='old',action='read')
1399 c      print *,"Processor",myrank," opened file IELEP"
1400       call getenv_loc('SIDEPAR',sidename)
1401       open (isidep,file=sidename,status='old',action='read')
1402 c      print *,"Processor",myrank," opened file ISIDEP"
1403 c      print *,"Processor",myrank," opened parameter files"
1404 #elif (defined G77)
1405       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1406       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1407 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1408 C Get parameter filenames and open the parameter files.
1409       call getenv_loc('BONDPAR',bondname)
1410       open (ibond,file=bondname,status='old')
1411       call getenv_loc('THETPAR',thetname)
1412       open (ithep,file=thetname,status='old')
1413 #ifndef CRYST_THETA
1414       call getenv_loc('THETPARPDB',thetname_pdb)
1415       open (ithep_pdb,file=thetname_pdb,status='old')
1416 #endif
1417       call getenv_loc('ROTPAR',rotname)
1418       open (irotam,file=rotname,status='old')
1419 #ifndef CRYST_SC
1420       call getenv_loc('ROTPARPDB',rotname_pdb)
1421       open (irotam_pdb,file=rotname_pdb,status='old')
1422 #endif
1423       call getenv_loc('TORPAR',torname)
1424       open (itorp,file=torname,status='old')
1425       call getenv_loc('TORDPAR',tordname)
1426       open (itordp,file=tordname,status='old')
1427       call getenv_loc('SCCORPAR',sccorname)
1428       open (isccor,file=sccorname,status='old')
1429       call getenv_loc('FOURIER',fouriername)
1430       open (ifourier,file=fouriername,status='old')
1431       call getenv_loc('ELEPAR',elename)
1432       open (ielep,file=elename,status='old')
1433       call getenv_loc('SIDEPAR',sidename)
1434       open (isidep,file=sidename,status='old')
1435 #else
1436       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1437      &  readonly)
1438        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1439 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1440 C Get parameter filenames and open the parameter files.
1441       call getenv_loc('BONDPAR',bondname)
1442       open (ibond,file=bondname,status='old',readonly)
1443       call getenv_loc('THETPAR',thetname)
1444       open (ithep,file=thetname,status='old',readonly)
1445 #ifndef CRYST_THETA
1446       call getenv_loc('THETPARPDB',thetname_pdb)
1447       print *,"thetname_pdb ",thetname_pdb
1448       open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1449       print *,ithep_pdb," opened"
1450 #endif
1451       call getenv_loc('ROTPAR',rotname)
1452       open (irotam,file=rotname,status='old',readonly)
1453 #ifndef CRYST_SC
1454       call getenv_loc('ROTPARPDB',rotname_pdb)
1455       open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1456 #endif
1457       call getenv_loc('TORPAR',torname)
1458       open (itorp,file=torname,status='old',readonly)
1459       call getenv_loc('TORDPAR',tordname)
1460       open (itordp,file=tordname,status='old',readonly)
1461       call getenv_loc('SCCORPAR',sccorname)
1462       open (isccor,file=sccorname,status='old',readonly)
1463       call getenv_loc('FOURIER',fouriername)
1464       open (ifourier,file=fouriername,status='old',readonly)
1465       call getenv_loc('ELEPAR',elename)
1466       open (ielep,file=elename,status='old',readonly)
1467       call getenv_loc('SIDEPAR',sidename)
1468       open (isidep,file=sidename,status='old',readonly)
1469 #endif
1470 #ifndef OLDSCP
1471 C
1472 C 8/9/01 In the newest version SCp interaction constants are read from a file
1473 C Use -DOLDSCP to use hard-coded constants instead.
1474 C
1475       call getenv_loc('SCPPAR',scpname)
1476 #if defined(WINIFL) || defined(WINPGI)
1477       open (iscpp,file=scpname,status='old',readonly,shared)
1478 #elif (defined CRAY)  || (defined AIX)
1479       open (iscpp,file=scpname,status='old',action='read')
1480 #elif (defined G77)
1481       open (iscpp,file=scpname,status='old')
1482 #else
1483       open (iscpp,file=scpname,status='old',readonly)
1484 #endif
1485 #endif
1486       call getenv_loc('PATTERN',patname)
1487 #if defined(WINIFL) || defined(WINPGI)
1488       open (icbase,file=patname,status='old',readonly,shared)
1489 #elif (defined CRAY)  || (defined AIX)
1490       open (icbase,file=patname,status='old',action='read')
1491 #elif (defined G77)
1492       open (icbase,file=patname,status='old')
1493 #else
1494       open (icbase,file=patname,status='old',readonly)
1495 #endif
1496 #ifdef MPI
1497 C Open output file only for CG processes
1498 c      print *,"Processor",myrank," fg_rank",fg_rank
1499       if (fg_rank.eq.0) then
1500
1501       if (nodes.eq.1) then
1502         npos=3
1503       else
1504         npos = dlog10(dfloat(nodes-1))+1
1505       endif
1506       if (npos.lt.3) npos=3
1507       write (liczba,'(i1)') npos
1508       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1509      &  //')'
1510       write (liczba,form) me
1511       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1512      &  liczba(:ilen(liczba))
1513       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1514      &  //'.int'
1515       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1516      &  //'.pdb'
1517       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1518      &  liczba(:ilen(liczba))//'.mol2'
1519       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1520      &  liczba(:ilen(liczba))//'.stat'
1521       if (lentmp.gt.0)
1522      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1523      &      //liczba(:ilen(liczba))//'.stat')
1524       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1525      &  //'.rst'
1526 c      if(usampl) then
1527 c          qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1528 c     & liczba(:ilen(liczba))//'.const'
1529 c      endif
1530
1531       endif
1532 #else
1533       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1534       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1535       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1536       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1537       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1538       if (lentmp.gt.0)
1539      & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1540      &    '.stat')
1541       rest2name=prefix(:ilen(prefix))//'.rst'
1542 c      if(usampl) then
1543 c         qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1544 c      endif
1545 #endif
1546 #if defined(AIX) || defined(PGI)
1547       if (me.eq.king .or. .not. out1file)
1548      &   open(iout,file=outname,status='unknown')
1549 #ifdef DEBUG
1550       if (fg_rank.gt.0) then
1551         write (liczba,'(i3.3)') myrank/nfgtasks
1552         write (ll,'(bz,i3.3)') fg_rank
1553         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1554      &   status='unknown')
1555       endif
1556 #endif
1557       if(me.eq.king) then
1558        open(igeom,file=intname,status='unknown',position='append')
1559        open(ipdb,file=pdbname,status='unknown')
1560        open(imol2,file=mol2name,status='unknown')
1561        open(istat,file=statname,status='unknown',position='append')
1562       else
1563 c1out       open(iout,file=outname,status='unknown')
1564       endif
1565 #else
1566       if (me.eq.king .or. .not.out1file)
1567      &    open(iout,file=outname,status='unknown')
1568 #ifdef DEBUG
1569       if (fg_rank.gt.0) then
1570         write (liczba,'(i3.3)') myrank/nfgtasks
1571         write (ll,'(bz,i3.3)') fg_rank
1572         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1573      &   status='unknown')
1574       endif
1575 #endif
1576       if(me.eq.king) then
1577        open(igeom,file=intname,status='unknown',access='append')
1578        open(ipdb,file=pdbname,status='unknown')
1579        open(imol2,file=mol2name,status='unknown')
1580        open(istat,file=statname,status='unknown',access='append')
1581       else
1582 c1out       open(iout,file=outname,status='unknown')
1583       endif
1584 #endif
1585       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1586       csa_seed=prefix(:lenpre)//'.CSA.seed'
1587       csa_history=prefix(:lenpre)//'.CSA.history'
1588       csa_bank=prefix(:lenpre)//'.CSA.bank'
1589       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1590       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1591       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1592 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1593       csa_int=prefix(:lenpre)//'.int'
1594       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1595       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1596       csa_in=prefix(:lenpre)//'.CSA.in'
1597 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1598 C Write file names
1599       if (me.eq.king)then
1600       write (iout,'(80(1h-))')
1601       write (iout,'(30x,a)') "FILE ASSIGNMENT"
1602       write (iout,'(80(1h-))')
1603       write (iout,*) "Input file                      : ",
1604      &  pref_orig(:ilen(pref_orig))//'.inp'
1605       write (iout,*) "Output file                     : ",
1606      &  outname(:ilen(outname))
1607       write (iout,*)
1608       write (iout,*) "Sidechain potential file        : ",
1609      &  sidename(:ilen(sidename))
1610 #ifndef OLDSCP
1611       write (iout,*) "SCp potential file              : ",
1612      &  scpname(:ilen(scpname))
1613 #endif
1614       write (iout,*) "Electrostatic potential file    : ",
1615      &  elename(:ilen(elename))
1616       write (iout,*) "Cumulant coefficient file       : ",
1617      &  fouriername(:ilen(fouriername))
1618       write (iout,*) "Torsional parameter file        : ",
1619      &  torname(:ilen(torname))
1620       write (iout,*) "Double torsional parameter file : ",
1621      &  tordname(:ilen(tordname))
1622       write (iout,*) "SCCOR parameter file : ",
1623      &  sccorname(:ilen(sccorname))
1624       write (iout,*) "Bond & inertia constant file    : ",
1625      &  bondname(:ilen(bondname))
1626       write (iout,*) "Bending parameter file          : ",
1627      &  thetname(:ilen(thetname))
1628       write (iout,*) "Rotamer parameter file          : ",
1629      &  rotname(:ilen(rotname))
1630       write (iout,*) "Threading database              : ",
1631      &  patname(:ilen(patname))
1632       if (lentmp.ne.0)
1633      &write (iout,*)" DIRTMP                          : ",
1634      &  tmpdir(:lentmp)
1635       write (iout,'(80(1h-))')
1636       endif
1637       return
1638       end
1639 c----------------------------------------------------------------------------
1640       subroutine card_concat(card)
1641       implicit real*8 (a-h,o-z)
1642       include 'DIMENSIONS'
1643       include 'COMMON.IOUNITS'
1644       character*(*) card
1645       character*80 karta,ucase
1646       external ilen
1647       read (inp,'(a)') karta
1648       karta=ucase(karta)
1649       card=' '
1650       do while (karta(80:80).eq.'&')
1651         card=card(:ilen(card)+1)//karta(:79)
1652         read (inp,'(a)') karta
1653         karta=ucase(karta)
1654       enddo
1655       card=card(:ilen(card)+1)//karta
1656       return
1657       end
1658
1659       subroutine read_dist_constr
1660       implicit real*8 (a-h,o-z)
1661       include 'DIMENSIONS'
1662 #ifdef MPI
1663       include 'mpif.h'
1664 #endif
1665       include 'COMMON.SETUP'
1666       include 'COMMON.CONTROL'
1667       include 'COMMON.CHAIN'
1668       include 'COMMON.IOUNITS'
1669       include 'COMMON.SBRIDGE'
1670       integer ifrag_(2,100),ipair_(2,100)
1671       double precision wfrag_(100),wpair_(100)
1672       character*500 controlcard
1673 c      write (iout,*) "Calling read_dist_constr"
1674 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1675 c      call flush(iout)
1676       call card_concat(controlcard)
1677       call readi(controlcard,"NFRAG",nfrag_,0)
1678       call readi(controlcard,"NPAIR",npair_,0)
1679       call readi(controlcard,"NDIST",ndist_,0)
1680       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1681       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1682       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1683       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1684       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1685 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1686 c      write (iout,*) "IFRAG"
1687 c      do i=1,nfrag_
1688 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1689 c      enddo
1690 c      write (iout,*) "IPAIR"
1691 c      do i=1,npair_
1692 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1693 c      enddo
1694       if (.not.refstr .and. nfrag.gt.0) then
1695         write (iout,*)
1696      &  "ERROR: no reference structure to compute distance restraints"
1697         write (iout,*)
1698      &  "Restraints must be specified explicitly (NDIST=number)"
1699         stop
1700       endif
1701       if (nfrag.lt.2 .and. npair.gt.0) then
1702         write (iout,*) "ERROR: Less than 2 fragments specified",
1703      &   " but distance restraints between pairs requested"
1704         stop
1705       endif
1706       call flush(iout)
1707       do i=1,nfrag_
1708         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1709         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1710      &    ifrag_(2,i)=nstart_sup+nsup-1
1711 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1712         call flush(iout)
1713         if (wfrag_(i).gt.0.0d0) then
1714         do j=ifrag_(1,i),ifrag_(2,i)-1
1715           do k=j+1,ifrag_(2,i)
1716 c            write (iout,*) "j",j," k",k
1717             ddjk=dist(j,k)
1718             if (constr_dist.eq.1) then
1719             nhpb=nhpb+1
1720             ihpb(nhpb)=j
1721             jhpb(nhpb)=k
1722               dhpb(nhpb)=ddjk
1723             forcon(nhpb)=wfrag_(i)
1724             else if (constr_dist.eq.2) then
1725               if (ddjk.le.dist_cut) then
1726                 nhpb=nhpb+1
1727                 ihpb(nhpb)=j
1728                 jhpb(nhpb)=k
1729                 dhpb(nhpb)=ddjk
1730                 forcon(nhpb)=wfrag_(i)
1731               endif
1732             else
1733               nhpb=nhpb+1
1734               ihpb(nhpb)=j
1735               jhpb(nhpb)=k
1736               dhpb(nhpb)=ddjk
1737               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1738             endif
1739 #ifdef MPI
1740             if (.not.out1file .or. me.eq.king)
1741      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1742      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1743 #else
1744             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1745      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1746 #endif
1747           enddo
1748         enddo
1749         endif
1750       enddo
1751       do i=1,npair_
1752         if (wpair_(i).gt.0.0d0) then
1753         ii = ipair_(1,i)
1754         jj = ipair_(2,i)
1755         if (ii.gt.jj) then
1756           itemp=ii
1757           ii=jj
1758           jj=itemp
1759         endif
1760         do j=ifrag_(1,ii),ifrag_(2,ii)
1761           do k=ifrag_(1,jj),ifrag_(2,jj)
1762             nhpb=nhpb+1
1763             ihpb(nhpb)=j
1764             jhpb(nhpb)=k
1765             forcon(nhpb)=wpair_(i)
1766             dhpb(nhpb)=dist(j,k)
1767 #ifdef MPI
1768             if (.not.out1file .or. me.eq.king)
1769      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1770      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1771 #else
1772             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1773      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1774 #endif
1775           enddo
1776         enddo
1777         endif
1778       enddo
1779       do i=1,ndist_
1780         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
1781      &     ibecarb(i),forcon(nhpb+1)
1782         if (forcon(nhpb+1).gt.0.0d0) then
1783           nhpb=nhpb+1
1784           if (ibecarb(i).gt.0) then
1785             ihpb(i)=ihpb(i)+nres
1786             jhpb(i)=jhpb(i)+nres
1787           endif
1788           if (dhpb(nhpb).eq.0.0d0)
1789      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1790         endif
1791       enddo
1792 #ifdef MPI
1793       if (.not.out1file .or. me.eq.king) then
1794 #endif
1795       do i=1,nhpb
1796           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
1797      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
1798       enddo
1799       call flush(iout)
1800 #ifdef MPI
1801       endif
1802 #endif
1803       return
1804       end
1805 c-------------------------------------------------------------------------------
1806 #ifdef WINIFL
1807       subroutine flush(iu)
1808       return
1809       end
1810 #endif
1811 #ifdef AIX
1812       subroutine flush(iu)
1813       call flush_(iu)
1814       return
1815       end
1816 #endif
1817 c------------------------------------------------------------------------------
1818       subroutine copy_to_tmp(source)
1819       include "DIMENSIONS"
1820       include "COMMON.IOUNITS"
1821       character*(*) source
1822       character* 256 tmpfile
1823       integer ilen
1824       external ilen
1825       logical ex
1826       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1827       inquire(file=tmpfile,exist=ex)
1828       if (ex) then
1829         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1830      &   " to temporary directory..."
1831         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1832         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1833       endif
1834       return
1835       end
1836 c------------------------------------------------------------------------------
1837       subroutine move_from_tmp(source)
1838       include "DIMENSIONS"
1839       include "COMMON.IOUNITS"
1840       character*(*) source
1841       integer ilen
1842       external ilen
1843       write (*,*) "Moving ",source(:ilen(source)),
1844      & " from temporary directory to working directory"
1845       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1846       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1847       return
1848       end
1849 c------------------------------------------------------------------------------
1850       subroutine random_init(seed)
1851 C
1852 C Initialize random number generator
1853 C
1854       implicit real*8 (a-h,o-z)
1855       include 'DIMENSIONS'
1856 #ifdef AMD64
1857       integer*8 iseedi8
1858 #endif
1859 #ifdef MPI
1860       include 'mpif.h'
1861       logical OKRandom, prng_restart
1862       real*8  r1
1863       integer iseed_array(4)
1864 #endif
1865       include 'COMMON.IOUNITS'
1866       include 'COMMON.TIME1'
1867 c      include 'COMMON.THREAD'
1868       include 'COMMON.SBRIDGE'
1869       include 'COMMON.CONTROL'
1870       include 'COMMON.MCM'
1871 c      include 'COMMON.MAP'
1872       include 'COMMON.HEADER'
1873 c      include 'COMMON.CSA'
1874       include 'COMMON.CHAIN'
1875 c      include 'COMMON.MUCA'
1876 c      include 'COMMON.MD'
1877       include 'COMMON.FFIELD'
1878       include 'COMMON.SETUP'
1879       iseed=-dint(dabs(seed))
1880       if (iseed.eq.0) then
1881         write (iout,'(/80(1h*)/20x,a/80(1h*))')
1882      &    'Random seed undefined. The program will stop.'
1883         write (*,'(/80(1h*)/20x,a/80(1h*))')
1884      &    'Random seed undefined. The program will stop.'
1885 #ifdef MPI
1886         call mpi_finalize(mpi_comm_world,ierr)
1887 #endif
1888         stop 'Bad random seed.'
1889       endif
1890 #ifdef MPI
1891       if (fg_rank.eq.0) then
1892       seed=seed*(me+1)+1
1893 #ifdef AMD64
1894       iseedi8=dint(seed)
1895       if(me.eq.king .or. .not. out1file)
1896      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1897       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1898       OKRandom = prng_restart(me,iseedi8)
1899 #else
1900       do i=1,4
1901        tmp=65536.0d0**(4-i)
1902        iseed_array(i) = dint(seed/tmp)
1903        seed=seed-iseed_array(i)*tmp
1904       enddo
1905       if(me.eq.king .or. .not. out1file)
1906      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1907      &                 (iseed_array(i),i=1,4)
1908       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1909      &                 (iseed_array(i),i=1,4)
1910       OKRandom = prng_restart(me,iseed_array)
1911 #endif
1912       if (OKRandom) then
1913         r1=ran_number(0.0D0,1.0D0)
1914         if(me.eq.king .or. .not. out1file)
1915      &   write (iout,*) 'ran_num',r1
1916         if (r1.lt.0.0d0) OKRandom=.false.
1917       endif
1918       if (.not.OKRandom) then
1919         write (iout,*) 'PRNG IS NOT WORKING!!!'
1920         print *,'PRNG IS NOT WORKING!!!'
1921         if (me.eq.0) then
1922          call flush(iout)
1923          call mpi_abort(mpi_comm_world,error_msg,ierr)
1924          stop
1925         else
1926          write (iout,*) 'too many processors for parallel prng'
1927          write (*,*) 'too many processors for parallel prng'
1928          call flush(iout)
1929          stop
1930         endif
1931       endif
1932       endif
1933 #else
1934       call vrndst(iseed)
1935       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
1936 #endif
1937       return
1938       end