2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
12 C Read force-field parameters except weights
14 C Read job setup parameters
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 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
20 C Read MD control parameters if reqjuired
21 if (modecalc.eq.12) call read_MDpar
22 C Read MREMD control parameters if required
23 if (modecalc.eq.14) then
27 C Read MUCA control parameters if required
28 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
35 cfmc if (modecalc.eq.10) call mcmfread
36 C Read molecule information, molecule geometry, energy-term weights, and
37 C restraints if requested
39 C Print restraint information
41 if (.not. out1file .or. me.eq.king) then
44 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
45 & " distance constraints have been imposed"
47 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
52 c print *,"Processor",myrank," leaves READRTNS"
55 C-------------------------------------------------------------------------------
56 subroutine read_control
60 implicit real*8 (a-h,o-z)
64 logical OKRandom, prng_restart
67 include 'COMMON.IOUNITS'
68 include 'COMMON.TIME1'
69 include 'COMMON.THREAD'
70 include 'COMMON.SBRIDGE'
71 include 'COMMON.CONTROL'
74 include 'COMMON.HEADER'
76 include 'COMMON.CHAIN'
79 include 'COMMON.FFIELD'
80 include 'COMMON.SETUP'
81 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
82 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
84 character*320 controlcard
89 read (INP,'(a)') titel
90 call card_concat(controlcard)
91 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
92 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
93 call reada(controlcard,'SEED',seed,0.0D0)
94 call random_init(seed)
95 C Set up the time limit (caution! The time must be input in minutes!)
96 read_cart=index(controlcard,'READ_CART').gt.0
97 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
98 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
99 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
100 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
101 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
102 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
103 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
104 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
105 call reada(controlcard,'DRMS',drms,0.1D0)
106 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
107 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
108 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
109 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
110 write (iout,'(a,f10.1)')'DRMS = ',drms
111 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
112 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
114 call readi(controlcard,'NZ_START',nz_start,0)
115 call readi(controlcard,'NZ_END',nz_end,0)
116 call readi(controlcard,'IZ_SC',iz_sc,0)
118 safety = 60.0d0*safety
121 call reada(controlcard,"T_BATH",t_bath,300.0d0)
122 minim=(index(controlcard,'MINIMIZE').gt.0)
123 dccart=(index(controlcard,'CART').gt.0)
124 overlapsc=(index(controlcard,'OVERLAP').gt.0)
125 overlapsc=.not.overlapsc
126 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
127 searchsc=.not.searchsc
128 sideadd=(index(controlcard,'SIDEADD').gt.0)
129 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
130 outpdb=(index(controlcard,'PDBOUT').gt.0)
131 outmol2=(index(controlcard,'MOL2OUT').gt.0)
132 pdbref=(index(controlcard,'PDBREF').gt.0)
133 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
134 indpdb=index(controlcard,'PDBSTART')
135 extconf=(index(controlcard,'EXTCONF').gt.0)
136 call readi(controlcard,'IPRINT',iprint,0)
137 call readi(controlcard,'MAXGEN',maxgen,10000)
138 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
139 call readi(controlcard,"KDIAG",kdiag,0)
140 call readi(controlcard,"RESCALE_MODE",rescale_mode,1)
141 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
142 & write (iout,*) "RESCALE_MODE",rescale_mode
143 split_ene=index(controlcard,'SPLIT_ENE').gt.0
144 if (index(controlcard,'REGULAR').gt.0.0D0) then
145 call reada(controlcard,'WEIDIS',weidis,0.1D0)
149 if (index(controlcard,'CHECKGRAD').gt.0) then
151 if (index(controlcard,'CART').gt.0) then
153 elseif (index(controlcard,'CARINT').gt.0) then
158 elseif (index(controlcard,'THREAD').gt.0) then
160 call readi(controlcard,'THREAD',nthread,0)
161 if (nthread.gt.0) then
162 call reada(controlcard,'WEIDIS',weidis,0.1D0)
165 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
166 stop 'Error termination in Read_Control.'
168 else if (index(controlcard,'MCMA').gt.0) then
170 else if (index(controlcard,'MCEE').gt.0) then
172 else if (index(controlcard,'MULTCONF').gt.0) then
174 else if (index(controlcard,'MAP').gt.0) then
176 call readi(controlcard,'MAP',nmap,0)
177 else if (index(controlcard,'CSA').gt.0) then
179 crc else if (index(controlcard,'ZSCORE').gt.0) then
181 crc ZSCORE is rm from UNRES, modecalc=9 is available
184 cfcm else if (index(controlcard,'MCMF').gt.0) then
186 else if (index(controlcard,'SOFTREG').gt.0) then
188 else if (index(controlcard,'CHECK_BOND').gt.0) then
190 else if (index(controlcard,'TEST').gt.0) then
192 else if (index(controlcard,'MD').gt.0) then
194 else if (index(controlcard,'RE ').gt.0) then
198 lmuca=index(controlcard,'MUCA').gt.0
199 call readi(controlcard,'MUCADYN',mucadyn,0)
200 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
201 if (lmuca .and. (me.eq.king .or. .not.out1file ))
203 write (iout,*) 'MUCADYN=',mucadyn
204 write (iout,*) 'MUCASMOOTH=',muca_smooth
207 iscode=index(controlcard,'ONE_LETTER')
208 indphi=index(controlcard,'PHI')
209 indback=index(controlcard,'BACK')
210 iranconf=index(controlcard,'RAND_CONF')
211 i2ndstr=index(controlcard,'USE_SEC_PRED')
212 gradout=index(controlcard,'GRADOUT').gt.0
213 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
215 if(me.eq.king.or..not.out1file)
216 & write (iout,'(2a)') diagmeth(kdiag),
217 & ' routine used to diagonalize matrices.'
220 c--------------------------------------------------------------------------
221 subroutine read_REMDpar
225 implicit real*8 (a-h,o-z)
227 include 'COMMON.IOUNITS'
228 include 'COMMON.TIME1'
231 include 'COMMON.LANGEVIN'
233 include 'COMMON.LANGEVIN.lang0'
235 include 'COMMON.INTERACT'
236 include 'COMMON.NAMES'
238 include 'COMMON.REMD'
239 include 'COMMON.CONTROL'
240 include 'COMMON.SETUP'
242 character*320 controlcard
243 character*3200 controlcard1
244 integer iremd_m_total
246 if(me.eq.king.or..not.out1file)
247 & write (iout,*) "REMD setup"
249 call card_concat(controlcard)
250 call readi(controlcard,"NREP",nrep,3)
251 call readi(controlcard,"NSTEX",nstex,1000)
252 call reada(controlcard,"RETMIN",retmin,10.0d0)
253 call reada(controlcard,"RETMAX",retmax,1000.0d0)
254 mremdsync=(index(controlcard,'SYNC').gt.0)
255 call readi(controlcard,"NSYN",i_sync_step,100)
256 restart1file=(index(controlcard,'REST1FILE').gt.0)
257 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
258 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
259 if(max_cache_traj_use.gt.max_cache_traj)
260 & max_cache_traj_use=max_cache_traj
261 if(me.eq.king.or..not.out1file) then
262 cd if (traj1file) then
263 crc caching is in testing - NTWX is not ignored
264 cd write (iout,*) "NTWX value is ignored"
265 cd write (iout,*) " trajectory is stored to one file by master"
266 cd write (iout,*) " before exchange at NSTEX intervals"
268 write (iout,*) "NREP= ",nrep
269 write (iout,*) "NSTEX= ",nstex
270 write (iout,*) "SYNC= ",mremdsync
271 write (iout,*) "NSYN= ",i_sync_step
272 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
275 t_exchange_only=(index(controlcard,'TONLY').gt.0)
276 call readi(controlcard,"HREMD",hremd,0)
277 if((me.eq.king.or..not.out1file).and.hremd.gt.0) then
278 write (iout,*) "Hamiltonian REMD with ",hremd," sets of weights"
280 if(usampl.and.hremd.gt.0) then
282 & "========== ERROR: USAMPL and HREMD cannot be used together"
284 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
291 if (index(controlcard,'TLIST').gt.0) then
293 call card_concat(controlcard1)
294 read(controlcard1,*) (remd_t(i),i=1,nrep)
295 if(me.eq.king.or..not.out1file)
296 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
299 if (index(controlcard,'MLIST').gt.0) then
301 call card_concat(controlcard1)
302 read(controlcard1,*) (remd_m(i),i=1,nrep)
303 if(me.eq.king.or..not.out1file) then
304 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
307 iremd_m_total=iremd_m_total+remd_m(i)
310 write (iout,*) 'Total number of replicas ',
311 & iremd_m_total*hremd
313 write (iout,*) 'Total number of replicas ',iremd_m_total
317 if(me.eq.king.or..not.out1file)
318 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
321 c--------------------------------------------------------------------------
322 subroutine read_MDpar
326 implicit real*8 (a-h,o-z)
328 include 'COMMON.IOUNITS'
329 include 'COMMON.TIME1'
332 include 'COMMON.LANGEVIN'
334 include 'COMMON.LANGEVIN.lang0'
336 include 'COMMON.INTERACT'
337 include 'COMMON.NAMES'
339 include 'COMMON.SETUP'
340 include 'COMMON.CONTROL'
341 include 'COMMON.SPLITELE'
343 character*320 controlcard
345 call card_concat(controlcard)
346 call readi(controlcard,"NSTEP",n_timestep,1000000)
347 call readi(controlcard,"NTWE",ntwe,100)
348 call readi(controlcard,"NTWX",ntwx,1000)
349 call reada(controlcard,"DT",d_time,1.0d-1)
350 call reada(controlcard,"DVMAX",dvmax,2.0d1)
351 call reada(controlcard,"DAMAX",damax,1.0d1)
352 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
353 call readi(controlcard,"LANG",lang,0)
354 RESPA = index(controlcard,"RESPA") .gt. 0
355 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
356 ntime_split0=ntime_split
357 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
358 ntime_split0=ntime_split
359 call reada(controlcard,"R_CUT",r_cut,2.0d0)
360 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
361 rest = index(controlcard,"REST").gt.0
362 tbf = index(controlcard,"TBF").gt.0
363 tnp = index(controlcard,"NOSEPOINCARE99").gt.0
364 tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0
365 tnh = index(controlcard,"NOSEHOOVER96").gt.0
366 if (RESPA.and.tnh)then
367 xiresp = index(controlcard,"XIRESP").gt.0
369 call reada(controlcard,"Q_NP",Q_np,0.1d0)
370 usampl = index(controlcard,"USAMPL").gt.0
372 mdpdb = index(controlcard,"MDPDB").gt.0
373 call reada(controlcard,"T_BATH",t_bath,300.0d0)
374 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
375 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
376 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
377 if (count_reset_moment.eq.0) count_reset_moment=1000000000
378 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
379 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
380 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
381 if (count_reset_vel.eq.0) count_reset_vel=1000000000
382 large = index(controlcard,"LARGE").gt.0
383 print_compon = index(controlcard,"PRINT_COMPON").gt.0
384 rattle = index(controlcard,"RATTLE").gt.0
385 c if performing umbrella sampling, fragments constrained are read from the fragment file
391 if(me.eq.king.or..not.out1file) then
393 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
395 write (iout,'(a)') "The units are:"
396 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
397 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
398 & " acceleration: angstrom/(48.9 fs)**2"
399 write (iout,'(a)') "energy: kcal/mol, temperature: K"
401 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
402 write (iout,'(a60,f10.5,a)')
403 & "Initial time step of numerical integration:",d_time,
405 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
407 write (iout,'(2a,i4,a)')
408 & "A-MTS algorithm used; initial time step for fast-varying",
409 & " short-range forces split into",ntime_split," steps."
410 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
411 & r_cut," lambda",rlamb
413 write (iout,'(2a,f10.5)')
414 & "Maximum acceleration threshold to reduce the time step",
415 & "/increase split number:",damax
416 write (iout,'(2a,f10.5)')
417 & "Maximum predicted energy drift to reduce the timestep",
418 & "/increase split number:",edriftmax
419 write (iout,'(a60,f10.5)')
420 & "Maximum velocity threshold to reduce velocities:",dvmax
421 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
422 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
423 if (rattle) write (iout,'(a60)')
424 & "Rattle algorithm used to constrain the virtual bonds"
428 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
429 call reada(controlcard,"RWAT",rwat,1.4d0)
430 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
431 surfarea=index(controlcard,"SURFAREA").gt.0
432 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
433 if(me.eq.king.or..not.out1file)then
434 write (iout,'(/a,$)') "Langevin dynamics calculation"
437 & " with direct integration of Langevin equations"
438 else if (lang.eq.2) then
439 write (iout,'(a/)') " with TINKER stochasic MD integrator"
440 else if (lang.eq.3) then
441 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
442 else if (lang.eq.4) then
443 write (iout,'(a/)') " in overdamped mode"
445 write (iout,'(//a,i5)')
446 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
449 write (iout,'(a60,f10.5)') "Temperature:",t_bath
450 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
451 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
452 write (iout,'(a60,f10.5)')
453 & "Scaling factor of the friction forces:",scal_fric
454 if (surfarea) write (iout,'(2a,i10,a)')
455 & "Friction coefficients will be scaled by solvent-accessible",
456 & " surface area every",reset_fricmat," steps."
458 c Calculate friction coefficients and bounds of stochastic forces
459 eta=6*pi*cPoise*etawat
460 if(me.eq.king.or..not.out1file)
461 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
463 gamp=scal_fric*(pstok+rwat)*eta
464 stdfp=dsqrt(2*Rb*t_bath/d_time)
466 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
467 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
469 if(me.eq.king.or..not.out1file)then
470 write (iout,'(/2a/)')
471 & "Radii of site types and friction coefficients and std's of",
472 & " stochastic forces of fully exposed sites"
473 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
475 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
476 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
480 if(me.eq.king.or..not.out1file)then
481 write (iout,'(a)') "Berendsen bath calculation"
482 write (iout,'(a60,f10.5)') "Temperature:",t_bath
483 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
485 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
486 & count_reset_moment," steps"
488 & write (iout,'(a,i10,a)')
489 & "Velocities will be reset at random every",count_reset_vel,
492 else if (tnp .or. tnp1 .or. tnh) then
493 if (tnp .or. tnp1) then
494 write (iout,'(a)') "Nose-Poincare bath calculation"
495 if (tnp) write (iout,'(a)')
496 & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird"
497 if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose"
499 write (iout,'(a)') "Nose-Hoover bath calculation"
500 write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al."
510 WDTI(i) = 1.0*d_time/nresn
517 write (iout,'(a)') "NVT-XI-RESPA algorithm"
519 write (iout,'(a)') "NVT-XO-RESPA algorithm"
522 WDTIi(i) = 1.0*d_time/nresn/ntime_split
530 write (iout,'(a60,f10.5)') "Temperature:",t_bath
531 write (iout,'(a60,f10.5)') "Q =",Q_np
533 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
534 & count_reset_moment," steps"
536 & write (iout,'(a,i10,a)')
537 & "Velocities will be reset at random every",count_reset_vel,
541 if(me.eq.king.or..not.out1file)
542 & write (iout,'(a31)') "Microcanonical mode calculation"
544 if(me.eq.king.or..not.out1file)then
545 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
547 write(iout,*) "MD running with constraints."
548 write(iout,*) "Equilibration time ", eq_time, " mtus."
549 write(iout,*) "Constraining ", nfrag," fragments."
550 write(iout,*) "Length of each fragment, weight and q0:"
552 write (iout,*) "Set of restraints #",iset
554 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
555 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
557 write(iout,*) "constraints between ", npair, "fragments."
558 write(iout,*) "constraint pairs, weights and q0:"
560 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
561 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
563 write(iout,*) "angle constraints within ", nfrag_back,
564 & "backbone fragments."
565 write(iout,*) "fragment, weights:"
567 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
568 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
569 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
572 iset=mod(kolor,nset)+1
575 if(me.eq.king.or..not.out1file)
576 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
579 c------------------------------------------------------------------------------
582 C Read molecular data.
584 implicit real*8 (a-h,o-z)
590 include 'COMMON.IOUNITS'
593 include 'COMMON.INTERACT'
594 include 'COMMON.LOCAL'
595 include 'COMMON.NAMES'
596 include 'COMMON.CHAIN'
597 include 'COMMON.FFIELD'
598 include 'COMMON.SBRIDGE'
599 include 'COMMON.HEADER'
600 include 'COMMON.CONTROL'
601 include 'COMMON.DBASE'
602 include 'COMMON.THREAD'
603 include 'COMMON.CONTACTS'
604 include 'COMMON.TORCNSTR'
605 include 'COMMON.TIME1'
606 include 'COMMON.BOUNDS'
608 include 'COMMON.REMD'
609 include 'COMMON.SETUP'
610 character*4 sequence(maxres)
612 double precision x(maxvar)
613 character*256 pdbfile
614 character*320 weightcard
615 character*80 weightcard_t,ucase
616 dimension itype_pdb(maxres)
617 common /pizda/ itype_pdb
618 logical seq_comp,fail
619 double precision energia(0:n_ene)
625 C Read weights of the subsequent energy terms.
638 if(me.eq.king.or..not.out1file) then
639 write (iout,*) 'Reading ',hremd,' sets of weights for HREMD'
640 write (iout,*) 'Current weights for processor ',
641 & me,' set ',i2set(me)
645 call card_concat(weightcard)
646 call reada(weightcard,'WLONG',wlong,1.0D0)
647 call reada(weightcard,'WSC',wsc,wlong)
648 call reada(weightcard,'WSCP',wscp,wlong)
649 call reada(weightcard,'WELEC',welec,1.0D0)
650 call reada(weightcard,'WVDWPP',wvdwpp,welec)
651 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
652 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
653 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
654 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
655 call reada(weightcard,'WTURN3',wturn3,1.0D0)
656 call reada(weightcard,'WTURN4',wturn4,1.0D0)
657 call reada(weightcard,'WTURN6',wturn6,1.0D0)
658 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
659 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
660 call reada(weightcard,'WBOND',wbond,1.0D0)
661 call reada(weightcard,'WTOR',wtor,1.0D0)
662 call reada(weightcard,'WTORD',wtor_d,1.0D0)
663 call reada(weightcard,'WANG',wang,1.0D0)
664 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
665 call reada(weightcard,'SCAL14',scal14,0.4D0)
666 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
667 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
668 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
669 call reada(weightcard,'TEMP0',temp0,300.0d0)
670 if (index(weightcard,'SOFT').gt.0) ipot=6
671 C 12/1/95 Added weight for the multi-body term WCORR
672 call reada(weightcard,'WCORRH',wcorr,1.0D0)
673 if (wcorr4.gt.0.0d0) wcorr=wcorr4
681 hweights(i,7)=wel_loc
684 hweights(i,10)=wturn6
686 hweights(i,12)=wscloc
688 hweights(i,14)=wtor_d
689 hweights(i,15)=wstrain
690 hweights(i,16)=wvdwpp
692 hweights(i,18)=scal14
693 hweights(i,21)=wsccor
698 weights(i)=hweights(i2set(me),i)
722 call card_concat(weightcard)
723 call reada(weightcard,'WLONG',wlong,1.0D0)
724 call reada(weightcard,'WSC',wsc,wlong)
725 call reada(weightcard,'WSCP',wscp,wlong)
726 call reada(weightcard,'WELEC',welec,1.0D0)
727 call reada(weightcard,'WVDWPP',wvdwpp,welec)
728 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
729 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
730 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
731 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
732 call reada(weightcard,'WTURN3',wturn3,1.0D0)
733 call reada(weightcard,'WTURN4',wturn4,1.0D0)
734 call reada(weightcard,'WTURN6',wturn6,1.0D0)
735 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
736 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
737 call reada(weightcard,'WBOND',wbond,1.0D0)
738 call reada(weightcard,'WTOR',wtor,1.0D0)
739 call reada(weightcard,'WTORD',wtor_d,1.0D0)
740 call reada(weightcard,'WANG',wang,1.0D0)
741 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
742 call reada(weightcard,'SCAL14',scal14,0.4D0)
743 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
744 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
745 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
746 call reada(weightcard,'TEMP0',temp0,300.0d0)
747 if (index(weightcard,'SOFT').gt.0) ipot=6
748 C 12/1/95 Added weight for the multi-body term WCORR
749 call reada(weightcard,'WCORRH',wcorr,1.0D0)
750 if (wcorr4.gt.0.0d0) wcorr=wcorr4
772 if(me.eq.king.or..not.out1file)
773 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
774 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
776 10 format (/'Energy-term weights (unscaled):'//
777 & 'WSCC= ',f10.6,' (SC-SC)'/
778 & 'WSCP= ',f10.6,' (SC-p)'/
779 & 'WELEC= ',f10.6,' (p-p electr)'/
780 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
781 & 'WBOND= ',f10.6,' (stretching)'/
782 & 'WANG= ',f10.6,' (bending)'/
783 & 'WSCLOC= ',f10.6,' (SC local)'/
784 & 'WTOR= ',f10.6,' (torsional)'/
785 & 'WTORD= ',f10.6,' (double torsional)'/
786 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
787 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
788 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
789 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
790 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
791 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
792 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
793 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
794 & 'WTURN6= ',f10.6,' (turns, 6th order)')
795 if(me.eq.king.or..not.out1file)then
796 if (wcorr4.gt.0.0d0) then
797 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
798 & 'between contact pairs of peptide groups'
799 write (iout,'(2(a,f5.3/))')
800 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
801 & 'Range of quenching the correlation terms:',2*delt_corr
802 else if (wcorr.gt.0.0d0) then
803 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
804 & 'between contact pairs of peptide groups'
806 write (iout,'(a,f8.3)')
807 & 'Scaling factor of 1,4 SC-p interactions:',scal14
808 write (iout,'(a,f8.3)')
809 & 'General scaling factor of SC-p interactions:',scalscp
811 r0_corr=cutoff_corr-delt_corr
813 aad(i,1)=scalscp*aad(i,1)
814 aad(i,2)=scalscp*aad(i,2)
815 bad(i,1)=scalscp*bad(i,1)
816 bad(i,2)=scalscp*bad(i,2)
818 call rescale_weights(t_bath)
819 if(me.eq.king.or..not.out1file)
820 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
821 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
823 22 format (/'Energy-term weights (scaled):'//
824 & 'WSCC= ',f10.6,' (SC-SC)'/
825 & 'WSCP= ',f10.6,' (SC-p)'/
826 & 'WELEC= ',f10.6,' (p-p electr)'/
827 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
828 & 'WBOND= ',f10.6,' (stretching)'/
829 & 'WANG= ',f10.6,' (bending)'/
830 & 'WSCLOC= ',f10.6,' (SC local)'/
831 & 'WTOR= ',f10.6,' (torsional)'/
832 & 'WTORD= ',f10.6,' (double torsional)'/
833 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
834 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
835 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
836 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
837 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
838 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
839 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
840 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
841 & 'WTURN6= ',f10.6,' (turns, 6th order)')
842 if(me.eq.king.or..not.out1file)
843 & write (iout,*) "Reference temperature for weights calculation:",
845 call reada(weightcard,"D0CM",d0cm,3.78d0)
846 call reada(weightcard,"AKCM",akcm,15.1d0)
847 call reada(weightcard,"AKTH",akth,11.0d0)
848 call reada(weightcard,"AKCT",akct,12.0d0)
849 call reada(weightcard,"V1SS",v1ss,-1.08d0)
850 call reada(weightcard,"V2SS",v2ss,7.61d0)
851 call reada(weightcard,"V3SS",v3ss,13.7d0)
852 call reada(weightcard,"EBR",ebr,-5.50D0)
853 if(me.eq.king.or..not.out1file) then
854 write (iout,*) "Parameters of the SS-bond potential:"
855 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
857 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
858 write (iout,*) "EBR",ebr
859 print *,'indpdb=',indpdb,' pdbref=',pdbref
861 if (indpdb.gt.0 .or. pdbref) then
862 read(inp,'(a)') pdbfile
863 if(me.eq.king.or..not.out1file)
864 & write (iout,'(2a)') 'PDB data will be read from file ',
865 & pdbfile(:ilen(pdbfile))
866 open(ipdbin,file=pdbfile,status='old',err=33)
868 33 write (iout,'(a)') 'Error opening PDB file.'
871 c print *,'Begin reading pdb data'
873 c print *,'Finished reading pdb data'
874 if(me.eq.king.or..not.out1file)
875 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
876 & ' nstart_sup=',nstart_sup
878 itype_pdb(i)=itype(i)
882 nct=nstart_sup+nsup-1
883 call contact(.false.,ncont_ref,icont_ref,co)
886 if(me.eq.king.or..not.out1file)
887 & write(iout,*)'Adding sidechains'
894 do while (fail.and.nsi.le.maxsi)
895 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
898 if(fail) write(iout,*)'Adding sidechain failed for res ',
899 & i,' after ',nsi,' trials'
904 if (indpdb.eq.0) then
905 C Read sequence if not taken from the pdb file.
907 c print *,'nres=',nres
908 if (iscode.gt.0) then
909 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
911 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
913 C Convert sequence to numeric code
915 itype(i)=rescode(i,sequence(i),iscode)
917 C Assign initial virtual bond lengths
923 vbld(i+nres)=dsc(itype(i))
924 vbld_inv(i+nres)=dsc_inv(itype(i))
925 c write (iout,*) "i",i," itype",itype(i),
926 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
930 c print '(20i4)',(itype(i),i=1,nres)
933 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
935 if (itype(i).eq.21) then
939 else if (itype(i+1).ne.20) then
941 else if (itype(i).ne.20) then
948 if(me.eq.king.or..not.out1file)then
949 write (iout,*) "ITEL"
951 write (iout,*) i,itype(i),itel(i)
953 print *,'Call Read_Bridge.'
956 C 8/13/98 Set limits to generating the dihedral angles
961 read (inp,*) ndih_constr
962 if (ndih_constr.gt.0) then
964 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
965 if(me.eq.king.or..not.out1file)then
967 & 'There are',ndih_constr,' constraints on phi angles.'
969 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
973 phi0(i)=deg2rad*phi0(i)
974 drange(i)=deg2rad*drange(i)
976 if(me.eq.king.or..not.out1file)
977 & write (iout,*) 'FTORS',ftors
980 phibound(1,ii) = phi0(i)-drange(i)
981 phibound(2,ii) = phi0(i)+drange(i)
988 write (iout,'(a)') 'Boundaries in phi angle sampling:'
990 write (iout,'(a3,i5,2f10.1)')
991 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
997 cd print *,'NNT=',NNT,' NCT=',NCT
998 if (itype(1).eq.21) nnt=2
999 if (itype(nres).eq.21) nct=nct-1
1001 if(me.eq.king.or..not.out1file)
1002 & write (iout,'(a,i3)') 'nsup=',nsup
1004 if (nsup.le.(nct-nnt+1)) then
1005 do i=0,nct-nnt+1-nsup
1006 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1012 & 'Error - sequences to be superposed do not match.'
1015 do i=0,nsup-(nct-nnt+1)
1016 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1018 nstart_sup=nstart_sup+i
1024 & 'Error - sequences to be superposed do not match.'
1027 if (nsup.eq.0) nsup=nct-nnt
1028 if (nstart_sup.eq.0) nstart_sup=nnt
1029 if (nstart_seq.eq.0) nstart_seq=nnt
1030 if(me.eq.king.or..not.out1file)
1031 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1032 & ' nstart_seq=',nstart_seq
1034 c--- Zscore rms -------
1035 if (nz_start.eq.0) nz_start=nnt
1036 if (nz_end.eq.0 .and. nsup.gt.0) then
1038 else if (nz_end.eq.0) then
1041 if(me.eq.king.or..not.out1file)then
1042 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1043 write (iout,*) 'IZ_SC=',iz_sc
1045 c----------------------
1048 if (.not.pdbref) then
1049 call read_angles(inp,*38)
1051 38 write (iout,'(a)') 'Error reading reference structure.'
1053 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1054 stop 'Error reading reference structure'
1058 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1067 call contact(.true.,ncont_ref,icont_ref,co)
1069 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1071 if (constr_dist.gt.0) call read_dist_constr
1072 c write (iout,*) "After read_dist_constr nhpb",nhpb
1074 if(me.eq.king.or..not.out1file)
1075 & write (iout,*) 'Contact order:',co
1077 if(me.eq.king.or..not.out1file)
1078 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1081 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1083 if(me.eq.king.or..not.out1file)
1084 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1085 & icont_ref(1,i),' ',
1086 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1090 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1091 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1092 & modecalc.ne.10) then
1093 C If input structure hasn't been supplied from the PDB file read or generate
1095 if (iranconf.eq.0 .and. .not. extconf) then
1096 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1097 & write (iout,'(a)') 'Initial geometry will be read in.'
1099 read(inp,'(8f10.5)',end=36,err=36)
1100 & ((c(l,k),l=1,3),k=1,nres),
1101 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1102 call int_from_cart1(.false.)
1105 dc(j,i)=c(j,i+1)-c(j,i)
1106 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1110 if (itype(i).ne.10) then
1112 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1113 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1119 call read_angles(inp,*36)
1122 36 write (iout,'(a)') 'Error reading angle file.'
1124 call mpi_finalize( MPI_COMM_WORLD,IERR )
1126 stop 'Error reading angle file.'
1128 else if (extconf) then
1129 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1130 & write (iout,'(a)') 'Extended chain initial geometry.'
1132 theta(i)=90d0*deg2rad
1135 phi(i)=180d0*deg2rad
1138 alph(i)=110d0*deg2rad
1141 omeg(i)=-120d0*deg2rad
1144 if(me.eq.king.or..not.out1file)
1145 & write (iout,'(a)') 'Random-generated initial geometry.'
1149 if (me.eq.king .or. fg_rank.eq.0 .and. (
1150 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1154 call gen_rand_conf(itmp,*30)
1156 30 write (iout,*) 'Failed to generate random conformation',
1157 & ', itrial=',itrial
1158 write (*,*) 'Processor:',me,
1159 & ' Failed to generate random conformation',
1169 write (iout,'(a,i3,a)') 'Processor:',me,
1170 & ' error in generating random conformation.'
1171 write (*,'(a,i3,a)') 'Processor:',me,
1172 & ' error in generating random conformation.'
1175 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1182 elseif (modecalc.eq.4) then
1183 read (inp,'(a)') intinname
1184 open (intin,file=intinname,status='old',err=333)
1185 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1186 & write (iout,'(a)') 'intinname',intinname
1187 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1189 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1191 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1193 stop 'Error opening angle file.'
1197 C Generate distance constraints, if the PDB structure is to be regularized.
1198 if (nthread.gt.0) then
1199 call read_threadbase
1202 if (me.eq.king .or. .not. out1file)
1204 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1205 write (iout,'(/a,i3,a)')
1206 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1207 write (iout,'(20i4)') (iss(i),i=1,ns)
1208 write (iout,'(/a/)') 'Pre-formed links are:'
1214 if (me.eq.king.or..not.out1file)
1215 & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1216 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1221 if (i2ndstr.gt.0) call secstrp2dihc
1222 c call geom_to_var(nvar,x)
1223 c call etotal(energia(0))
1224 c call enerprint(energia(0))
1225 c call briefout(0,etot)
1227 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1228 cd write (iout,'(a)') 'Variable list:'
1229 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1231 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1232 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1233 & 'Processor',myrank,': end reading molecular data.'
1237 c--------------------------------------------------------------------------
1238 logical function seq_comp(itypea,itypeb,length)
1240 integer length,itypea(length),itypeb(length)
1243 if (itypea(i).ne.itypeb(i)) then
1251 c-----------------------------------------------------------------------------
1252 subroutine read_bridge
1253 C Read information about disulfide bridges.
1254 implicit real*8 (a-h,o-z)
1255 include 'DIMENSIONS'
1259 include 'COMMON.IOUNITS'
1260 include 'COMMON.GEO'
1261 include 'COMMON.VAR'
1262 include 'COMMON.INTERACT'
1263 include 'COMMON.LOCAL'
1264 include 'COMMON.NAMES'
1265 include 'COMMON.CHAIN'
1266 include 'COMMON.FFIELD'
1267 include 'COMMON.SBRIDGE'
1268 include 'COMMON.HEADER'
1269 include 'COMMON.CONTROL'
1270 include 'COMMON.DBASE'
1271 include 'COMMON.THREAD'
1272 include 'COMMON.TIME1'
1273 include 'COMMON.SETUP'
1274 C Read bridging residues.
1275 read (inp,*) ns,(iss(i),i=1,ns)
1277 if(me.eq.king.or..not.out1file)
1278 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1279 C Check whether the specified bridging residues are cystines.
1281 if (itype(iss(i)).ne.1) then
1282 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1283 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
1284 & ' can form a disulfide bridge?!!!'
1285 write (*,'(2a,i3,a)')
1286 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
1287 & ' can form a disulfide bridge?!!!'
1289 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1294 C Read preformed bridges.
1296 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1297 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1300 C Check if the residues involved in bridges are in the specified list of
1301 C bridging residues.
1304 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1305 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1306 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1307 & ' contains residues present in other pairs.'
1308 write (*,'(a,i3,a)') 'Disulfide pair',i,
1309 & ' contains residues present in other pairs.'
1311 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1317 if (ihpb(i).eq.iss(j)) goto 10
1319 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1322 if (jhpb(i).eq.iss(j)) goto 20
1324 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1330 ihpb(i)=ihpb(i)+nres
1331 jhpb(i)=jhpb(i)+nres
1337 c----------------------------------------------------------------------------
1338 subroutine read_x(kanal,*)
1339 implicit real*8 (a-h,o-z)
1340 include 'DIMENSIONS'
1341 include 'COMMON.GEO'
1342 include 'COMMON.VAR'
1343 include 'COMMON.CHAIN'
1344 include 'COMMON.IOUNITS'
1345 include 'COMMON.CONTROL'
1346 include 'COMMON.LOCAL'
1347 include 'COMMON.INTERACT'
1348 c Read coordinates from input
1350 read(kanal,'(8f10.5)',end=10,err=10)
1351 & ((c(l,k),l=1,3),k=1,nres),
1352 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1355 c(j,2*nres)=c(j,nres)
1357 call int_from_cart1(.false.)
1360 dc(j,i)=c(j,i+1)-c(j,i)
1361 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1365 if (itype(i).ne.10) then
1367 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1368 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1376 c----------------------------------------------------------------------------
1377 subroutine read_threadbase
1378 implicit real*8 (a-h,o-z)
1379 include 'DIMENSIONS'
1380 include 'COMMON.IOUNITS'
1381 include 'COMMON.GEO'
1382 include 'COMMON.VAR'
1383 include 'COMMON.INTERACT'
1384 include 'COMMON.LOCAL'
1385 include 'COMMON.NAMES'
1386 include 'COMMON.CHAIN'
1387 include 'COMMON.FFIELD'
1388 include 'COMMON.SBRIDGE'
1389 include 'COMMON.HEADER'
1390 include 'COMMON.CONTROL'
1391 include 'COMMON.DBASE'
1392 include 'COMMON.THREAD'
1393 include 'COMMON.TIME1'
1394 C Read pattern database for threading.
1395 read (icbase,*) nseq
1397 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1398 & nres_base(2,i),nres_base(3,i)
1399 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1401 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1402 c & nres_base(2,i),nres_base(3,i)
1403 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1407 if (weidis.eq.0.0D0) weidis=0.1D0
1416 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1417 write (iout,'(a,i5)') 'nexcl: ',nexcl
1418 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1421 c------------------------------------------------------------------------------
1422 subroutine setup_var
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 Set up variable list.
1445 if (itype(i).ne.10) then
1447 ialph(i,1)=nvar+nside
1451 if (indphi.gt.0) then
1453 else if (indback.gt.0) then
1458 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1461 c----------------------------------------------------------------------------
1462 subroutine gen_dist_constr
1463 C Generate CA distance constraints.
1464 implicit real*8 (a-h,o-z)
1465 include 'DIMENSIONS'
1466 include 'COMMON.IOUNITS'
1467 include 'COMMON.GEO'
1468 include 'COMMON.VAR'
1469 include 'COMMON.INTERACT'
1470 include 'COMMON.LOCAL'
1471 include 'COMMON.NAMES'
1472 include 'COMMON.CHAIN'
1473 include 'COMMON.FFIELD'
1474 include 'COMMON.SBRIDGE'
1475 include 'COMMON.HEADER'
1476 include 'COMMON.CONTROL'
1477 include 'COMMON.DBASE'
1478 include 'COMMON.THREAD'
1479 include 'COMMON.TIME1'
1480 dimension itype_pdb(maxres)
1481 common /pizda/ itype_pdb
1483 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1484 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1485 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1487 do i=nstart_sup,nstart_sup+nsup-1
1488 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1489 cd & ' seq_pdb', restyp(itype_pdb(i))
1490 do j=i+2,nstart_sup+nsup-1
1492 ihpb(nhpb)=i+nstart_seq-nstart_sup
1493 jhpb(nhpb)=j+nstart_seq-nstart_sup
1495 dhpb(nhpb)=dist(i,j)
1498 cd write (iout,'(a)') 'Distance constraints:'
1503 cd if (ii.gt.nres) then
1508 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1509 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1510 cd & dhpb(i),forcon(i)
1514 c----------------------------------------------------------------------------
1516 implicit real*8 (a-h,o-z)
1517 include 'DIMENSIONS'
1518 include 'COMMON.MAP'
1519 include 'COMMON.IOUNITS'
1520 character*3 angid(4) /'THE','PHI','ALP','OME'/
1521 character*80 mapcard,ucase
1523 read (inp,'(a)') mapcard
1524 mapcard=ucase(mapcard)
1525 if (index(mapcard,'PHI').gt.0) then
1527 else if (index(mapcard,'THE').gt.0) then
1529 else if (index(mapcard,'ALP').gt.0) then
1531 else if (index(mapcard,'OME').gt.0) then
1534 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1535 stop 'Error - illegal variable spec in MAP card.'
1537 call readi (mapcard,'RES1',res1(imap),0)
1538 call readi (mapcard,'RES2',res2(imap),0)
1539 if (res1(imap).eq.0) then
1540 res1(imap)=res2(imap)
1541 else if (res2(imap).eq.0) then
1542 res2(imap)=res1(imap)
1544 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1546 & 'Error - illegal definition of variable group in MAP.'
1547 stop 'Error - illegal definition of variable group in MAP.'
1549 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1550 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1551 call readi(mapcard,'NSTEP',nstep(imap),0)
1552 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1554 & 'Illegal boundary and/or step size specification in MAP.'
1555 stop 'Illegal boundary and/or step size specification in MAP.'
1560 c----------------------------------------------------------------------------
1562 implicit real*8 (a-h,o-z)
1563 include 'DIMENSIONS'
1564 include 'COMMON.IOUNITS'
1565 include 'COMMON.GEO'
1566 include 'COMMON.CSA'
1567 include 'COMMON.BANK'
1568 include 'COMMON.CONTROL'
1570 character*620 mcmcard
1571 call card_concat(mcmcard)
1573 call readi(mcmcard,'NCONF',nconf,50)
1574 call readi(mcmcard,'NADD',nadd,0)
1575 call readi(mcmcard,'JSTART',jstart,1)
1576 call readi(mcmcard,'JEND',jend,1)
1577 call readi(mcmcard,'NSTMAX',nstmax,500000)
1578 call readi(mcmcard,'N0',n0,1)
1579 call readi(mcmcard,'N1',n1,6)
1580 call readi(mcmcard,'N2',n2,4)
1581 call readi(mcmcard,'N3',n3,0)
1582 call readi(mcmcard,'N4',n4,0)
1583 call readi(mcmcard,'N5',n5,0)
1584 call readi(mcmcard,'N6',n6,10)
1585 call readi(mcmcard,'N7',n7,0)
1586 call readi(mcmcard,'N8',n8,0)
1587 call readi(mcmcard,'N9',n9,0)
1588 call readi(mcmcard,'N14',n14,0)
1589 call readi(mcmcard,'N15',n15,0)
1590 call readi(mcmcard,'N16',n16,0)
1591 call readi(mcmcard,'N17',n17,0)
1592 call readi(mcmcard,'N18',n18,0)
1594 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1596 call readi(mcmcard,'NDIFF',ndiff,2)
1597 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1598 call readi(mcmcard,'IS1',is1,1)
1599 call readi(mcmcard,'IS2',is2,8)
1600 call readi(mcmcard,'NRAN0',nran0,4)
1601 call readi(mcmcard,'NRAN1',nran1,2)
1602 call readi(mcmcard,'IRR',irr,1)
1603 call readi(mcmcard,'NSEED',nseed,20)
1604 call readi(mcmcard,'NTOTAL',ntotal,10000)
1605 call reada(mcmcard,'CUT1',cut1,2.0d0)
1606 call reada(mcmcard,'CUT2',cut2,5.0d0)
1607 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1608 call readi(mcmcard,'ICMAX',icmax,3)
1609 call readi(mcmcard,'IRESTART',irestart,0)
1610 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1613 call reada(mcmcard,'DELE',dele,20.0d0)
1614 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1615 call readi(mcmcard,'IREF',iref,0)
1616 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1617 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1618 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1619 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1620 write (iout,*) "NCONF_IN",nconf_in
1623 c----------------------------------------------------------------------------
1624 cfmc subroutine mcmfread
1625 cfmc implicit real*8 (a-h,o-z)
1626 cfmc include 'DIMENSIONS'
1627 cfmc include 'COMMON.MCMF'
1628 cfmc include 'COMMON.IOUNITS'
1629 cfmc include 'COMMON.GEO'
1630 cfmc character*80 ucase
1631 cfmc character*620 mcmcard
1632 cfmc call card_concat(mcmcard)
1634 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1635 cfmc write(iout,*)'MAXRANT=',maxrant
1636 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1637 cfmc write(iout,*)'MAXFAM=',maxfam
1638 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1639 cfmc write(iout,*)'NNET1=',nnet1
1640 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1641 cfmc write(iout,*)'NNET2=',nnet2
1642 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1643 cfmc write(iout,*)'NNET3=',nnet3
1644 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1645 cfmc write(iout,*)'ILASTT=',ilastt
1646 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1647 cfmc write(iout,*)'MAXSTR=',maxstr
1648 cfmc maxstr_f=maxstr/maxfam
1649 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1650 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1651 cfmc write(iout,*)'NMCMF=',nmcmf
1652 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1653 cfmc write(iout,*)'IFOCUS=',ifocus
1654 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1655 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1656 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1657 cfmc write(iout,*)'INTPRT=',intprt
1658 cfmc call readi(mcmcard,'IPRT',iprt,100)
1659 cfmc write(iout,*)'IPRT=',iprt
1660 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1661 cfmc write(iout,*)'IMAXTR=',imaxtr
1662 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1663 cfmc write(iout,*)'MAXEVEN=',maxeven
1664 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1665 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1666 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1667 cfmc write(iout,*)'INIMIN=',inimin
1668 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1669 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1670 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1671 cfmc write(iout,*)'NTHREAD=',nthread
1672 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1673 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1674 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1675 cfmc write(iout,*)'MAXPERT=',maxpert
1676 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1677 cfmc write(iout,*)'IRMSD=',irmsd
1678 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1679 cfmc write(iout,*)'DENEMIN=',denemin
1680 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1681 cfmc write(iout,*)'RCUT1S=',rcut1s
1682 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1683 cfmc write(iout,*)'RCUT1E=',rcut1e
1684 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1685 cfmc write(iout,*)'RCUT2S=',rcut2s
1686 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1687 cfmc write(iout,*)'RCUT2E=',rcut2e
1688 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1689 cfmc write(iout,*)'DPERT1=',d_pert1
1690 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1691 cfmc write(iout,*)'DPERT1A=',d_pert1a
1692 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1693 cfmc write(iout,*)'DPERT2=',d_pert2
1694 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1695 cfmc write(iout,*)'DPERT2A=',d_pert2a
1696 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1697 cfmc write(iout,*)'DPERT2B=',d_pert2b
1698 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1699 cfmc write(iout,*)'DPERT2C=',d_pert2c
1700 cfmc d_pert1=deg2rad*d_pert1
1701 cfmc d_pert1a=deg2rad*d_pert1a
1702 cfmc d_pert2=deg2rad*d_pert2
1703 cfmc d_pert2a=deg2rad*d_pert2a
1704 cfmc d_pert2b=deg2rad*d_pert2b
1705 cfmc d_pert2c=deg2rad*d_pert2c
1706 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1707 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1708 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1709 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1710 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1711 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1712 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1713 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1714 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1715 cfmc write(iout,*)'RCUTINI=',rcutini
1716 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1717 cfmc write(iout,*)'GRAT=',grat
1718 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1719 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1723 c----------------------------------------------------------------------------
1725 implicit real*8 (a-h,o-z)
1726 include 'DIMENSIONS'
1727 include 'COMMON.MCM'
1728 include 'COMMON.MCE'
1729 include 'COMMON.IOUNITS'
1731 character*320 mcmcard
1732 call card_concat(mcmcard)
1733 call readi(mcmcard,'MAXACC',maxacc,100)
1734 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1735 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1736 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1737 call readi(mcmcard,'MAXREPM',maxrepm,200)
1738 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1739 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1740 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1741 call reada(mcmcard,'E_UP',e_up,5.0D0)
1742 call reada(mcmcard,'DELTE',delte,0.1D0)
1743 call readi(mcmcard,'NSWEEP',nsweep,5)
1744 call readi(mcmcard,'NSTEPH',nsteph,0)
1745 call readi(mcmcard,'NSTEPC',nstepc,0)
1746 call reada(mcmcard,'TMIN',tmin,298.0D0)
1747 call reada(mcmcard,'TMAX',tmax,298.0D0)
1748 call readi(mcmcard,'NWINDOW',nwindow,0)
1749 call readi(mcmcard,'PRINT_MC',print_mc,0)
1750 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1751 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1752 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1753 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1754 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1755 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1756 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1757 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1758 if (nwindow.gt.0) then
1759 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1761 winlen(i)=winend(i)-winstart(i)+1
1764 if (tmax.lt.tmin) tmax=tmin
1765 if (tmax.eq.tmin) then
1769 if (nstepc.gt.0 .and. nsteph.gt.0) then
1770 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1771 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1773 C Probabilities of different move types
1774 sumpro_type(0)=0.0D0
1775 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1776 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1777 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1778 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1779 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1780 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1781 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1783 print *,'i',i,' sumprotype',sumpro_type(i)
1784 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1785 print *,'i',i,' sumprotype',sumpro_type(i)
1789 c----------------------------------------------------------------------------
1790 subroutine read_minim
1791 implicit real*8 (a-h,o-z)
1792 include 'DIMENSIONS'
1793 include 'COMMON.MINIM'
1794 include 'COMMON.IOUNITS'
1796 character*320 minimcard
1797 call card_concat(minimcard)
1798 call readi(minimcard,'MAXMIN',maxmin,2000)
1799 call readi(minimcard,'MAXFUN',maxfun,5000)
1800 call readi(minimcard,'MINMIN',minmin,maxmin)
1801 call readi(minimcard,'MINFUN',minfun,maxmin)
1802 call reada(minimcard,'TOLF',tolf,1.0D-2)
1803 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1804 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1805 & 'Options in energy minimization:'
1806 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1807 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1808 & 'MinMin:',MinMin,' MinFun:',MinFun,
1809 & ' TolF:',TolF,' RTolF:',RTolF
1812 c----------------------------------------------------------------------------
1813 subroutine read_angles(kanal,*)
1814 implicit real*8 (a-h,o-z)
1815 include 'DIMENSIONS'
1816 include 'COMMON.GEO'
1817 include 'COMMON.VAR'
1818 include 'COMMON.CHAIN'
1819 include 'COMMON.IOUNITS'
1820 include 'COMMON.CONTROL'
1821 c Read angles from input
1823 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1824 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1825 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1826 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1829 c 9/7/01 avoid 180 deg valence angle
1830 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1832 theta(i)=deg2rad*theta(i)
1833 phi(i)=deg2rad*phi(i)
1834 alph(i)=deg2rad*alph(i)
1835 omeg(i)=deg2rad*omeg(i)
1840 c----------------------------------------------------------------------------
1841 subroutine reada(rekord,lancuch,wartosc,default)
1843 character*(*) rekord,lancuch
1844 double precision wartosc,default
1847 iread=index(rekord,lancuch)
1848 if (iread.eq.0) then
1852 iread=iread+ilen(lancuch)+1
1853 read (rekord(iread:),*,err=10,end=10) wartosc
1858 c----------------------------------------------------------------------------
1859 subroutine readi(rekord,lancuch,wartosc,default)
1861 character*(*) rekord,lancuch
1862 integer wartosc,default
1865 iread=index(rekord,lancuch)
1866 if (iread.eq.0) then
1870 iread=iread+ilen(lancuch)+1
1871 read (rekord(iread:),*,err=10,end=10) wartosc
1876 c----------------------------------------------------------------------------
1877 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1880 integer tablica(dim),default
1881 character*(*) rekord,lancuch
1888 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1889 if (iread.eq.0) return
1890 iread=iread+ilen(lancuch)+1
1891 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1894 c----------------------------------------------------------------------------
1895 subroutine multreada(rekord,lancuch,tablica,dim,default)
1898 double precision tablica(dim),default
1899 character*(*) rekord,lancuch
1906 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1907 if (iread.eq.0) return
1908 iread=iread+ilen(lancuch)+1
1909 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1912 c----------------------------------------------------------------------------
1913 subroutine openunits
1914 implicit real*8 (a-h,o-z)
1915 include 'DIMENSIONS'
1918 character*16 form,nodename
1921 include 'COMMON.SETUP'
1922 include 'COMMON.IOUNITS'
1924 include 'COMMON.CONTROL'
1925 integer lenpre,lenpot,ilen,lentmp
1927 character*3 out1file_text,ucase
1930 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1931 call getenv_loc("PREFIX",prefix)
1933 call getenv_loc("POT",pot)
1934 call getenv_loc("DIRTMP",tmpdir)
1935 call getenv_loc("CURDIR",curdir)
1936 call getenv_loc("OUT1FILE",out1file_text)
1937 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1938 out1file_text=ucase(out1file_text)
1939 if (out1file_text(1:1).eq."Y") then
1942 out1file=fg_rank.gt.0
1947 if (lentmp.gt.0) then
1948 write (*,'(80(1h!))')
1949 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1950 write (*,'(80(1h!))')
1951 write (*,*)"All output files will be on node /tmp directory."
1953 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1954 if (me.eq.king) then
1955 write (*,*) "The master node is ",nodename
1956 else if (fg_rank.eq.0) then
1957 write (*,*) "I am the CG slave node ",nodename
1959 write (*,*) "I am the FG slave node ",nodename
1962 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1963 lenpre = lentmp+lenpre+1
1965 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1966 C Get the names and open the input files
1967 #if defined(WINIFL) || defined(WINPGI)
1968 open(1,file=pref_orig(:ilen(pref_orig))//
1969 & '.inp',status='old',readonly,shared)
1970 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1971 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1972 C Get parameter filenames and open the parameter files.
1973 call getenv_loc('BONDPAR',bondname)
1974 open (ibond,file=bondname,status='old',readonly,shared)
1975 call getenv_loc('THETPAR',thetname)
1976 open (ithep,file=thetname,status='old',readonly,shared)
1977 call getenv_loc('ROTPAR',rotname)
1978 open (irotam,file=rotname,status='old',readonly,shared)
1979 call getenv_loc('TORPAR',torname)
1980 open (itorp,file=torname,status='old',readonly,shared)
1981 call getenv_loc('TORDPAR',tordname)
1982 open (itordp,file=tordname,status='old',readonly,shared)
1983 call getenv_loc('FOURIER',fouriername)
1984 open (ifourier,file=fouriername,status='old',readonly,shared)
1985 call getenv_loc('ELEPAR',elename)
1986 open (ielep,file=elename,status='old',readonly,shared)
1987 call getenv_loc('SIDEPAR',sidename)
1988 open (isidep,file=sidename,status='old',readonly,shared)
1989 #elif (defined CRAY) || (defined AIX)
1990 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1992 c print *,"Processor",myrank," opened file 1"
1993 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1994 c print *,"Processor",myrank," opened file 9"
1995 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1996 C Get parameter filenames and open the parameter files.
1997 call getenv_loc('BONDPAR',bondname)
1998 open (ibond,file=bondname,status='old',action='read')
1999 c print *,"Processor",myrank," opened file IBOND"
2000 call getenv_loc('THETPAR',thetname)
2001 open (ithep,file=thetname,status='old',action='read')
2002 c print *,"Processor",myrank," opened file ITHEP"
2003 call getenv_loc('ROTPAR',rotname)
2004 open (irotam,file=rotname,status='old',action='read')
2005 c print *,"Processor",myrank," opened file IROTAM"
2006 call getenv_loc('TORPAR',torname)
2007 open (itorp,file=torname,status='old',action='read')
2008 c print *,"Processor",myrank," opened file ITORP"
2009 call getenv_loc('TORDPAR',tordname)
2010 open (itordp,file=tordname,status='old',action='read')
2011 c print *,"Processor",myrank," opened file ITORDP"
2012 call getenv_loc('SCCORPAR',sccorname)
2013 open (isccor,file=sccorname,status='old',action='read')
2014 c print *,"Processor",myrank," opened file ISCCOR"
2015 call getenv_loc('FOURIER',fouriername)
2016 open (ifourier,file=fouriername,status='old',action='read')
2017 c print *,"Processor",myrank," opened file IFOURIER"
2018 call getenv_loc('ELEPAR',elename)
2019 open (ielep,file=elename,status='old',action='read')
2020 c print *,"Processor",myrank," opened file IELEP"
2021 call getenv_loc('SIDEPAR',sidename)
2022 open (isidep,file=sidename,status='old',action='read')
2023 c print *,"Processor",myrank," opened file ISIDEP"
2024 c print *,"Processor",myrank," opened parameter files"
2026 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2027 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2028 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2029 C Get parameter filenames and open the parameter files.
2030 call getenv_loc('BONDPAR',bondname)
2031 open (ibond,file=bondname,status='old')
2032 call getenv_loc('THETPAR',thetname)
2033 open (ithep,file=thetname,status='old')
2034 call getenv_loc('ROTPAR',rotname)
2035 open (irotam,file=rotname,status='old')
2036 call getenv_loc('TORPAR',torname)
2037 open (itorp,file=torname,status='old')
2038 call getenv_loc('TORDPAR',tordname)
2039 open (itordp,file=tordname,status='old')
2040 call getenv_loc('SCCORPAR',sccorname)
2041 open (isccor,file=sccorname,status='old')
2042 call getenv_loc('FOURIER',fouriername)
2043 open (ifourier,file=fouriername,status='old')
2044 call getenv_loc('ELEPAR',elename)
2045 open (ielep,file=elename,status='old')
2046 call getenv_loc('SIDEPAR',sidename)
2047 open (isidep,file=sidename,status='old')
2049 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2051 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2052 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2053 C Get parameter filenames and open the parameter files.
2054 call getenv_loc('BONDPAR',bondname)
2055 open (ibond,file=bondname,status='old',readonly)
2056 call getenv_loc('THETPAR',thetname)
2057 open (ithep,file=thetname,status='old',readonly)
2058 call getenv_loc('ROTPAR',rotname)
2059 open (irotam,file=rotname,status='old',readonly)
2060 call getenv_loc('TORPAR',torname)
2061 open (itorp,file=torname,status='old',readonly)
2062 call getenv_loc('TORDPAR',tordname)
2063 open (itordp,file=tordname,status='old',readonly)
2064 call getenv_loc('SCCORPAR',sccorname)
2065 open (isccor,file=sccorname,status='old',readonly)
2066 call getenv_loc('FOURIER',fouriername)
2067 open (ifourier,file=fouriername,status='old',readonly)
2068 call getenv_loc('ELEPAR',elename)
2069 open (ielep,file=elename,status='old',readonly)
2070 call getenv_loc('SIDEPAR',sidename)
2071 open (isidep,file=sidename,status='old',readonly)
2075 C 8/9/01 In the newest version SCp interaction constants are read from a file
2076 C Use -DOLDSCP to use hard-coded constants instead.
2078 call getenv_loc('SCPPAR',scpname)
2079 #if defined(WINIFL) || defined(WINPGI)
2080 open (iscpp,file=scpname,status='old',readonly,shared)
2081 #elif (defined CRAY) || (defined AIX)
2082 open (iscpp,file=scpname,status='old',action='read')
2084 open (iscpp,file=scpname,status='old')
2086 open (iscpp,file=scpname,status='old',readonly)
2089 call getenv_loc('PATTERN',patname)
2090 #if defined(WINIFL) || defined(WINPGI)
2091 open (icbase,file=patname,status='old',readonly,shared)
2092 #elif (defined CRAY) || (defined AIX)
2093 open (icbase,file=patname,status='old',action='read')
2095 open (icbase,file=patname,status='old')
2097 open (icbase,file=patname,status='old',readonly)
2100 C Open output file only for CG processes
2101 c print *,"Processor",myrank," fg_rank",fg_rank
2102 if (fg_rank.eq.0) then
2104 if (nodes.eq.1) then
2107 npos = dlog10(dfloat(nodes-1))+1
2109 if (npos.lt.3) npos=3
2110 write (liczba,'(i1)') npos
2111 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2113 write (liczba,form) me
2114 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2115 & liczba(:ilen(liczba))
2116 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2118 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2120 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2121 & liczba(:ilen(liczba))//'.mol2'
2122 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2123 & liczba(:ilen(liczba))//'.stat'
2125 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2126 & //liczba(:ilen(liczba))//'.stat')
2127 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2130 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2131 & liczba(:ilen(liczba))//'.const'
2136 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2137 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2138 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2139 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2140 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2142 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2144 rest2name=prefix(:ilen(prefix))//'.rst'
2146 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2149 #if defined(AIX) || defined(PGI)
2150 if (me.eq.king .or. .not. out1file)
2151 & open(iout,file=outname,status='unknown')
2153 if (fg_rank.gt.0) then
2154 write (liczba,'(i3.3)') myrank/nfgtasks
2155 write (ll,'(bz,i3.3)') fg_rank
2156 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2161 open(igeom,file=intname,status='unknown',position='append')
2162 open(ipdb,file=pdbname,status='unknown')
2163 open(imol2,file=mol2name,status='unknown')
2164 open(istat,file=statname,status='unknown',position='append')
2166 c1out open(iout,file=outname,status='unknown')
2169 if (me.eq.king .or. .not.out1file)
2170 & open(iout,file=outname,status='unknown')
2172 if (fg_rank.gt.0) then
2173 write (liczba,'(i3.3)') myrank/nfgtasks
2174 write (ll,'(bz,i3.3)') fg_rank
2175 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2180 open(igeom,file=intname,status='unknown',access='append')
2181 open(ipdb,file=pdbname,status='unknown')
2182 open(imol2,file=mol2name,status='unknown')
2183 open(istat,file=statname,status='unknown',access='append')
2185 c1out open(iout,file=outname,status='unknown')
2188 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2189 csa_seed=prefix(:lenpre)//'.CSA.seed'
2190 csa_history=prefix(:lenpre)//'.CSA.history'
2191 csa_bank=prefix(:lenpre)//'.CSA.bank'
2192 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2193 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2194 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2195 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2196 csa_int=prefix(:lenpre)//'.int'
2197 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2198 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2199 csa_in=prefix(:lenpre)//'.CSA.in'
2200 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2203 write (iout,'(80(1h-))')
2204 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2205 write (iout,'(80(1h-))')
2206 write (iout,*) "Input file : ",
2207 & pref_orig(:ilen(pref_orig))//'.inp'
2208 write (iout,*) "Output file : ",
2209 & outname(:ilen(outname))
2211 write (iout,*) "Sidechain potential file : ",
2212 & sidename(:ilen(sidename))
2214 write (iout,*) "SCp potential file : ",
2215 & scpname(:ilen(scpname))
2217 write (iout,*) "Electrostatic potential file : ",
2218 & elename(:ilen(elename))
2219 write (iout,*) "Cumulant coefficient file : ",
2220 & fouriername(:ilen(fouriername))
2221 write (iout,*) "Torsional parameter file : ",
2222 & torname(:ilen(torname))
2223 write (iout,*) "Double torsional parameter file : ",
2224 & tordname(:ilen(tordname))
2225 write (iout,*) "SCCOR parameter file : ",
2226 & sccorname(:ilen(sccorname))
2227 write (iout,*) "Bond & inertia constant file : ",
2228 & bondname(:ilen(bondname))
2229 write (iout,*) "Bending parameter file : ",
2230 & thetname(:ilen(thetname))
2231 write (iout,*) "Rotamer parameter file : ",
2232 & rotname(:ilen(rotname))
2233 write (iout,*) "Threading database : ",
2234 & patname(:ilen(patname))
2236 &write (iout,*)" DIRTMP : ",
2238 write (iout,'(80(1h-))')
2242 c----------------------------------------------------------------------------
2243 subroutine card_concat(card)
2244 implicit real*8 (a-h,o-z)
2245 include 'DIMENSIONS'
2246 include 'COMMON.IOUNITS'
2248 character*80 karta,ucase
2250 read (inp,'(a)') karta
2253 do while (karta(80:80).eq.'&')
2254 card=card(:ilen(card)+1)//karta(:79)
2255 read (inp,'(a)') karta
2258 card=card(:ilen(card)+1)//karta
2261 c----------------------------------------------------------------------------------
2263 implicit real*8 (a-h,o-z)
2264 include 'DIMENSIONS'
2265 include 'COMMON.CHAIN'
2266 include 'COMMON.IOUNITS'
2268 open(irest2,file=rest2name,status='unknown')
2269 read(irest2,*) totT,EK,potE,totE,t_bath
2271 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2274 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2277 read (irest2,*) iset
2282 c---------------------------------------------------------------------------------
2283 subroutine read_fragments
2284 implicit real*8 (a-h,o-z)
2285 include 'DIMENSIONS'
2289 include 'COMMON.SETUP'
2290 include 'COMMON.CHAIN'
2291 include 'COMMON.IOUNITS'
2293 include 'COMMON.CONTROL'
2294 read(inp,*) nset,nfrag,npair,nfrag_back
2295 if(me.eq.king.or..not.out1file)
2296 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2297 & " nfrag_back",nfrag_back
2299 read(inp,*) mset(iset)
2301 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2303 if(me.eq.king.or..not.out1file)
2304 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2305 & ifrag(2,i,iset), qinfrag(i,iset)
2308 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2310 if(me.eq.king.or..not.out1file)
2311 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2312 & ipair(2,i,iset), qinpair(i,iset)
2315 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2316 & wfrag_back(3,i,iset),
2317 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2318 if(me.eq.king.or..not.out1file)
2319 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2320 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2325 c-------------------------------------------------------------------------------
2326 subroutine read_dist_constr
2327 implicit real*8 (a-h,o-z)
2328 include 'DIMENSIONS'
2332 include 'COMMON.SETUP'
2333 include 'COMMON.CONTROL'
2334 include 'COMMON.CHAIN'
2335 include 'COMMON.IOUNITS'
2336 include 'COMMON.SBRIDGE'
2337 integer ifrag_(2,100),ipair_(2,100)
2338 double precision wfrag_(100),wpair_(100)
2339 character*500 controlcard
2340 c write (iout,*) "Calling read_dist_constr"
2341 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2343 call card_concat(controlcard)
2344 call readi(controlcard,"NFRAG",nfrag_,0)
2345 call readi(controlcard,"NPAIR",npair_,0)
2346 call readi(controlcard,"NDIST",ndist_,0)
2347 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2348 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2349 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2350 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2351 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2352 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2353 c write (iout,*) "IFRAG"
2355 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2357 c write (iout,*) "IPAIR"
2359 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2363 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2364 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2365 & ifrag_(2,i)=nstart_sup+nsup-1
2366 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2368 if (wfrag_(i).gt.0.0d0) then
2369 do j=ifrag_(1,i),ifrag_(2,i)-1
2370 do k=j+1,ifrag_(2,i)
2371 write (iout,*) "j",j," k",k
2373 if (constr_dist.eq.1) then
2378 forcon(nhpb)=wfrag_(i)
2379 else if (constr_dist.eq.2) then
2380 if (ddjk.le.dist_cut) then
2385 forcon(nhpb)=wfrag_(i)
2392 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2395 if (.not.out1file .or. me.eq.king)
2396 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2397 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2399 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2400 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2407 if (wpair_(i).gt.0.0d0) then
2415 do j=ifrag_(1,ii),ifrag_(2,ii)
2416 do k=ifrag_(1,jj),ifrag_(2,jj)
2420 forcon(nhpb)=wpair_(i)
2421 dhpb(nhpb)=dist(j,k)
2423 if (.not.out1file .or. me.eq.king)
2424 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2425 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2427 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2428 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2435 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2436 if (forcon(nhpb+1).gt.0.0d0) then
2438 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2440 if (.not.out1file .or. me.eq.king)
2441 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2442 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2444 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2445 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2452 c-------------------------------------------------------------------------------
2454 subroutine flush(iu)
2459 subroutine flush(iu)
2464 c------------------------------------------------------------------------------
2465 subroutine copy_to_tmp(source)
2466 include "DIMENSIONS"
2467 include "COMMON.IOUNITS"
2468 character*(*) source
2469 character* 256 tmpfile
2473 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2474 inquire(file=tmpfile,exist=ex)
2476 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2477 & " to temporary directory..."
2478 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2479 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2483 c------------------------------------------------------------------------------
2484 subroutine move_from_tmp(source)
2485 include "DIMENSIONS"
2486 include "COMMON.IOUNITS"
2487 character*(*) source
2490 write (*,*) "Moving ",source(:ilen(source)),
2491 & " from temporary directory to working directory"
2492 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2493 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2496 c------------------------------------------------------------------------------
2497 subroutine random_init(seed)
2499 C Initialize random number generator
2501 implicit real*8 (a-h,o-z)
2502 include 'DIMENSIONS'
2508 logical OKRandom, prng_restart
2510 integer iseed_array(4)
2512 include 'COMMON.IOUNITS'
2513 include 'COMMON.TIME1'
2514 include 'COMMON.THREAD'
2515 include 'COMMON.SBRIDGE'
2516 include 'COMMON.CONTROL'
2517 include 'COMMON.MCM'
2518 include 'COMMON.MAP'
2519 include 'COMMON.HEADER'
2520 include 'COMMON.CSA'
2521 include 'COMMON.CHAIN'
2522 include 'COMMON.MUCA'
2524 include 'COMMON.FFIELD'
2525 include 'COMMON.SETUP'
2526 iseed=-dint(dabs(seed))
2527 if (iseed.eq.0) then
2528 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2529 & 'Random seed undefined. The program will stop.'
2530 write (*,'(/80(1h*)/20x,a/80(1h*))')
2531 & 'Random seed undefined. The program will stop.'
2533 call mpi_finalize(mpi_comm_world,ierr)
2535 stop 'Bad random seed.'
2538 if (fg_rank.eq.0) then
2542 if(me.eq.king .or. .not. out1file)
2543 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2544 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2545 OKRandom = prng_restart(me,iseedi8)
2548 tmp=65536.0d0**(4-i)
2549 iseed_array(i) = dint(seed/tmp)
2550 seed=seed-iseed_array(i)*tmp
2552 if(me.eq.king .or. .not. out1file)
2553 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2554 & (iseed_array(i),i=1,4)
2555 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2556 & (iseed_array(i),i=1,4)
2557 OKRandom = prng_restart(me,iseed_array)
2560 r1=ran_number(0.0D0,1.0D0)
2561 if(me.eq.king .or. .not. out1file)
2562 & write (iout,*) 'ran_num',r1
2563 if (r1.lt.0.0d0) OKRandom=.false.
2565 if (.not.OKRandom) then
2566 write (iout,*) 'PRNG IS NOT WORKING!!!'
2567 print *,'PRNG IS NOT WORKING!!!'
2570 call mpi_abort(mpi_comm_world,error_msg,ierr)
2573 write (iout,*) 'too many processors for parallel prng'
2574 write (*,*) 'too many processors for parallel prng'
2582 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)