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.INTERACT'
81 include 'COMMON.SETUP'
82 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
83 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
85 character*320 controlcard
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 readi(controlcard,'SYM',symetr,1)
100 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
101 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
102 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
103 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
104 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
105 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
106 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
107 call reada(controlcard,'DRMS',drms,0.1D0)
108 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
109 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
110 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
111 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
112 write (iout,'(a,f10.1)')'DRMS = ',drms
113 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
114 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
116 call readi(controlcard,'NZ_START',nz_start,0)
117 call readi(controlcard,'NZ_END',nz_end,0)
118 call readi(controlcard,'IZ_SC',iz_sc,0)
120 safety = 60.0d0*safety
123 call reada(controlcard,"T_BATH",t_bath,300.0d0)
124 minim=(index(controlcard,'MINIMIZE').gt.0)
125 dccart=(index(controlcard,'CART').gt.0)
126 overlapsc=(index(controlcard,'OVERLAP').gt.0)
127 overlapsc=.not.overlapsc
128 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
129 searchsc=.not.searchsc
130 sideadd=(index(controlcard,'SIDEADD').gt.0)
131 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
132 outpdb=(index(controlcard,'PDBOUT').gt.0)
133 outmol2=(index(controlcard,'MOL2OUT').gt.0)
134 pdbref=(index(controlcard,'PDBREF').gt.0)
135 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
136 indpdb=index(controlcard,'PDBSTART')
137 extconf=(index(controlcard,'EXTCONF').gt.0)
138 call readi(controlcard,'IPRINT',iprint,0)
139 call readi(controlcard,'MAXGEN',maxgen,10000)
140 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
141 call readi(controlcard,"KDIAG",kdiag,0)
142 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
143 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
144 & write (iout,*) "RESCALE_MODE",rescale_mode
145 split_ene=index(controlcard,'SPLIT_ENE').gt.0
146 if (index(controlcard,'REGULAR').gt.0.0D0) then
147 call reada(controlcard,'WEIDIS',weidis,0.1D0)
151 if (index(controlcard,'CHECKGRAD').gt.0) then
153 if (index(controlcard,'CART').gt.0) then
155 elseif (index(controlcard,'CARINT').gt.0) then
160 elseif (index(controlcard,'THREAD').gt.0) then
162 call readi(controlcard,'THREAD',nthread,0)
163 if (nthread.gt.0) then
164 call reada(controlcard,'WEIDIS',weidis,0.1D0)
167 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
168 stop 'Error termination in Read_Control.'
170 else if (index(controlcard,'MCMA').gt.0) then
172 else if (index(controlcard,'MCEE').gt.0) then
174 else if (index(controlcard,'MULTCONF').gt.0) then
176 else if (index(controlcard,'MAP').gt.0) then
178 call readi(controlcard,'MAP',nmap,0)
179 else if (index(controlcard,'CSA').gt.0) then
181 crc else if (index(controlcard,'ZSCORE').gt.0) then
183 crc ZSCORE is rm from UNRES, modecalc=9 is available
186 cfcm else if (index(controlcard,'MCMF').gt.0) then
188 else if (index(controlcard,'SOFTREG').gt.0) then
190 else if (index(controlcard,'CHECK_BOND').gt.0) then
192 else if (index(controlcard,'TEST').gt.0) then
194 else if (index(controlcard,'MD').gt.0) then
196 else if (index(controlcard,'RE ').gt.0) then
200 lmuca=index(controlcard,'MUCA').gt.0
201 call readi(controlcard,'MUCADYN',mucadyn,0)
202 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
203 if (lmuca .and. (me.eq.king .or. .not.out1file ))
205 write (iout,*) 'MUCADYN=',mucadyn
206 write (iout,*) 'MUCASMOOTH=',muca_smooth
209 iscode=index(controlcard,'ONE_LETTER')
210 indphi=index(controlcard,'PHI')
211 indback=index(controlcard,'BACK')
212 iranconf=index(controlcard,'RAND_CONF')
213 i2ndstr=index(controlcard,'USE_SEC_PRED')
214 gradout=index(controlcard,'GRADOUT').gt.0
215 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
216 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
217 if (me.eq.king .or. .not.out1file )
218 & write (iout,*) "DISTCHAINMAX",distchainmax
220 if(me.eq.king.or..not.out1file)
221 & write (iout,'(2a)') diagmeth(kdiag),
222 & ' routine used to diagonalize matrices.'
225 c--------------------------------------------------------------------------
226 subroutine read_REMDpar
230 implicit real*8 (a-h,o-z)
232 include 'COMMON.IOUNITS'
233 include 'COMMON.TIME1'
236 include 'COMMON.LANGEVIN'
238 include 'COMMON.LANGEVIN.lang0'
240 include 'COMMON.INTERACT'
241 include 'COMMON.NAMES'
243 include 'COMMON.REMD'
244 include 'COMMON.CONTROL'
245 include 'COMMON.SETUP'
247 character*320 controlcard
248 character*3200 controlcard1
249 integer iremd_m_total
251 if(me.eq.king.or..not.out1file)
252 & write (iout,*) "REMD setup"
254 call card_concat(controlcard)
255 call readi(controlcard,"NREP",nrep,3)
256 call readi(controlcard,"NSTEX",nstex,1000)
257 call reada(controlcard,"RETMIN",retmin,10.0d0)
258 call reada(controlcard,"RETMAX",retmax,1000.0d0)
259 mremdsync=(index(controlcard,'SYNC').gt.0)
260 call readi(controlcard,"NSYN",i_sync_step,100)
261 restart1file=(index(controlcard,'REST1FILE').gt.0)
262 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
263 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
264 if(max_cache_traj_use.gt.max_cache_traj)
265 & max_cache_traj_use=max_cache_traj
266 if(me.eq.king.or..not.out1file) then
267 cd if (traj1file) then
268 crc caching is in testing - NTWX is not ignored
269 cd write (iout,*) "NTWX value is ignored"
270 cd write (iout,*) " trajectory is stored to one file by master"
271 cd write (iout,*) " before exchange at NSTEX intervals"
273 write (iout,*) "NREP= ",nrep
274 write (iout,*) "NSTEX= ",nstex
275 write (iout,*) "SYNC= ",mremdsync
276 write (iout,*) "NSYN= ",i_sync_step
277 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
280 if (index(controlcard,'TLIST').gt.0) then
282 call card_concat(controlcard1)
283 read(controlcard1,*) (remd_t(i),i=1,nrep)
284 if(me.eq.king.or..not.out1file)
285 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
288 if (index(controlcard,'MLIST').gt.0) then
290 call card_concat(controlcard1)
291 read(controlcard1,*) (remd_m(i),i=1,nrep)
292 if(me.eq.king.or..not.out1file) then
293 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
296 iremd_m_total=iremd_m_total+remd_m(i)
298 write (iout,*) 'Total number of replicas ',iremd_m_total
301 if(me.eq.king.or..not.out1file)
302 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
305 c--------------------------------------------------------------------------
306 subroutine read_MDpar
310 implicit real*8 (a-h,o-z)
312 include 'COMMON.IOUNITS'
313 include 'COMMON.TIME1'
316 include 'COMMON.LANGEVIN'
318 include 'COMMON.LANGEVIN.lang0'
320 include 'COMMON.INTERACT'
321 include 'COMMON.NAMES'
323 include 'COMMON.SETUP'
324 include 'COMMON.CONTROL'
325 include 'COMMON.SPLITELE'
327 character*320 controlcard
329 call card_concat(controlcard)
330 call readi(controlcard,"NSTEP",n_timestep,1000000)
331 call readi(controlcard,"NTWE",ntwe,100)
332 call readi(controlcard,"NTWX",ntwx,1000)
333 call reada(controlcard,"DT",d_time,1.0d-1)
334 call reada(controlcard,"DVMAX",dvmax,2.0d1)
335 call reada(controlcard,"DAMAX",damax,1.0d1)
336 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
337 call readi(controlcard,"LANG",lang,0)
338 RESPA = index(controlcard,"RESPA") .gt. 0
339 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
340 ntime_split0=ntime_split
341 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
342 ntime_split0=ntime_split
343 call reada(controlcard,"R_CUT",r_cut,2.0d0)
344 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
345 rest = index(controlcard,"REST").gt.0
346 tbf = index(controlcard,"TBF").gt.0
347 usampl = index(controlcard,"USAMPL").gt.0
349 mdpdb = index(controlcard,"MDPDB").gt.0
350 call reada(controlcard,"T_BATH",t_bath,300.0d0)
351 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
352 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
353 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
354 if (count_reset_moment.eq.0) count_reset_moment=1000000000
355 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
356 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
357 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
358 if (count_reset_vel.eq.0) count_reset_vel=1000000000
359 large = index(controlcard,"LARGE").gt.0
360 print_compon = index(controlcard,"PRINT_COMPON").gt.0
361 rattle = index(controlcard,"RATTLE").gt.0
362 c if performing umbrella sampling, fragments constrained are read from the fragment file
368 if(me.eq.king.or..not.out1file) then
370 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
372 write (iout,'(a)') "The units are:"
373 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
374 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
375 & " acceleration: angstrom/(48.9 fs)**2"
376 write (iout,'(a)') "energy: kcal/mol, temperature: K"
378 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
379 write (iout,'(a60,f10.5,a)')
380 & "Initial time step of numerical integration:",d_time,
382 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
384 write (iout,'(2a,i4,a)')
385 & "A-MTS algorithm used; initial time step for fast-varying",
386 & " short-range forces split into",ntime_split," steps."
387 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
388 & r_cut," lambda",rlamb
390 write (iout,'(2a,f10.5)')
391 & "Maximum acceleration threshold to reduce the time step",
392 & "/increase split number:",damax
393 write (iout,'(2a,f10.5)')
394 & "Maximum predicted energy drift to reduce the timestep",
395 & "/increase split number:",edriftmax
396 write (iout,'(a60,f10.5)')
397 & "Maximum velocity threshold to reduce velocities:",dvmax
398 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
399 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
400 if (rattle) write (iout,'(a60)')
401 & "Rattle algorithm used to constrain the virtual bonds"
405 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
406 call reada(controlcard,"RWAT",rwat,1.4d0)
407 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
408 surfarea=index(controlcard,"SURFAREA").gt.0
409 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
410 if(me.eq.king.or..not.out1file)then
411 write (iout,'(/a,$)') "Langevin dynamics calculation"
414 & " with direct integration of Langevin equations"
415 else if (lang.eq.2) then
416 write (iout,'(a/)') " with TINKER stochasic MD integrator"
417 else if (lang.eq.3) then
418 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
419 else if (lang.eq.4) then
420 write (iout,'(a/)') " in overdamped mode"
422 write (iout,'(//a,i5)')
423 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
426 write (iout,'(a60,f10.5)') "Temperature:",t_bath
427 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
428 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
429 write (iout,'(a60,f10.5)')
430 & "Scaling factor of the friction forces:",scal_fric
431 if (surfarea) write (iout,'(2a,i10,a)')
432 & "Friction coefficients will be scaled by solvent-accessible",
433 & " surface area every",reset_fricmat," steps."
435 c Calculate friction coefficients and bounds of stochastic forces
436 eta=6*pi*cPoise*etawat
437 if(me.eq.king.or..not.out1file)
438 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
440 gamp=scal_fric*(pstok+rwat)*eta
441 stdfp=dsqrt(2*Rb*t_bath/d_time)
443 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
444 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
446 if(me.eq.king.or..not.out1file)then
447 write (iout,'(/2a/)')
448 & "Radii of site types and friction coefficients and std's of",
449 & " stochastic forces of fully exposed sites"
450 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
452 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
453 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
457 if(me.eq.king.or..not.out1file)then
458 write (iout,'(a)') "Berendsen bath calculation"
459 write (iout,'(a60,f10.5)') "Temperature:",t_bath
460 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
462 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
463 & count_reset_moment," steps"
465 & write (iout,'(a,i10,a)')
466 & "Velocities will be reset at random every",count_reset_vel,
470 if(me.eq.king.or..not.out1file)
471 & write (iout,'(a31)') "Microcanonical mode calculation"
473 if(me.eq.king.or..not.out1file)then
474 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
476 write(iout,*) "MD running with constraints."
477 write(iout,*) "Equilibration time ", eq_time, " mtus."
478 write(iout,*) "Constraining ", nfrag," fragments."
479 write(iout,*) "Length of each fragment, weight and q0:"
481 write (iout,*) "Set of restraints #",iset
483 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
484 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
486 write(iout,*) "constraints between ", npair, "fragments."
487 write(iout,*) "constraint pairs, weights and q0:"
489 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
490 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
492 write(iout,*) "angle constraints within ", nfrag_back,
493 & "backbone fragments."
494 write(iout,*) "fragment, weights:"
496 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
497 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
498 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
501 iset=mod(kolor,nset)+1
504 if(me.eq.king.or..not.out1file)
505 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
508 c------------------------------------------------------------------------------
511 C Read molecular data.
513 implicit real*8 (a-h,o-z)
519 include 'COMMON.IOUNITS'
522 include 'COMMON.INTERACT'
523 include 'COMMON.LOCAL'
524 include 'COMMON.NAMES'
525 include 'COMMON.CHAIN'
526 include 'COMMON.FFIELD'
527 include 'COMMON.SBRIDGE'
528 include 'COMMON.HEADER'
529 include 'COMMON.CONTROL'
530 include 'COMMON.DBASE'
531 include 'COMMON.THREAD'
532 include 'COMMON.CONTACTS'
533 include 'COMMON.TORCNSTR'
534 include 'COMMON.TIME1'
535 include 'COMMON.BOUNDS'
537 include 'COMMON.SETUP'
538 character*4 sequence(maxres)
540 double precision x(maxvar)
541 character*256 pdbfile
542 character*400 weightcard
543 character*80 weightcard_t,ucase
544 dimension itype_pdb(maxres)
545 common /pizda/ itype_pdb
546 logical seq_comp,fail
547 double precision energia(0:n_ene)
553 C Read weights of the subsequent energy terms.
554 call card_concat(weightcard)
555 call reada(weightcard,'WLONG',wlong,1.0D0)
556 call reada(weightcard,'WSC',wsc,wlong)
557 call reada(weightcard,'WSCP',wscp,wlong)
558 call reada(weightcard,'WELEC',welec,1.0D0)
559 call reada(weightcard,'WVDWPP',wvdwpp,welec)
560 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
561 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
562 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
563 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
564 call reada(weightcard,'WTURN3',wturn3,1.0D0)
565 call reada(weightcard,'WTURN4',wturn4,1.0D0)
566 call reada(weightcard,'WTURN6',wturn6,1.0D0)
567 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
568 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
569 call reada(weightcard,'WBOND',wbond,1.0D0)
570 call reada(weightcard,'WTOR',wtor,1.0D0)
571 call reada(weightcard,'WTORD',wtor_d,1.0D0)
572 call reada(weightcard,'WANG',wang,1.0D0)
573 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
574 call reada(weightcard,'SCAL14',scal14,0.4D0)
575 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
576 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
577 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
578 call reada(weightcard,'TEMP0',temp0,300.0d0)
579 if (index(weightcard,'SOFT').gt.0) ipot=6
580 C 12/1/95 Added weight for the multi-body term WCORR
581 call reada(weightcard,'WCORRH',wcorr,1.0D0)
582 if (wcorr4.gt.0.0d0) wcorr=wcorr4
602 if(me.eq.king.or..not.out1file)
603 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
604 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
606 10 format (/'Energy-term weights (unscaled):'//
607 & 'WSCC= ',f10.6,' (SC-SC)'/
608 & 'WSCP= ',f10.6,' (SC-p)'/
609 & 'WELEC= ',f10.6,' (p-p electr)'/
610 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
611 & 'WBOND= ',f10.6,' (stretching)'/
612 & 'WANG= ',f10.6,' (bending)'/
613 & 'WSCLOC= ',f10.6,' (SC local)'/
614 & 'WTOR= ',f10.6,' (torsional)'/
615 & 'WTORD= ',f10.6,' (double torsional)'/
616 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
617 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
618 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
619 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
620 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
621 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
622 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
623 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
624 & 'WTURN6= ',f10.6,' (turns, 6th order)')
625 if(me.eq.king.or..not.out1file)then
626 if (wcorr4.gt.0.0d0) then
627 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
628 & 'between contact pairs of peptide groups'
629 write (iout,'(2(a,f5.3/))')
630 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
631 & 'Range of quenching the correlation terms:',2*delt_corr
632 else if (wcorr.gt.0.0d0) then
633 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
634 & 'between contact pairs of peptide groups'
636 write (iout,'(a,f8.3)')
637 & 'Scaling factor of 1,4 SC-p interactions:',scal14
638 write (iout,'(a,f8.3)')
639 & 'General scaling factor of SC-p interactions:',scalscp
641 r0_corr=cutoff_corr-delt_corr
643 aad(i,1)=scalscp*aad(i,1)
644 aad(i,2)=scalscp*aad(i,2)
645 bad(i,1)=scalscp*bad(i,1)
646 bad(i,2)=scalscp*bad(i,2)
648 call rescale_weights(t_bath)
649 if(me.eq.king.or..not.out1file)
650 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
651 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
653 22 format (/'Energy-term weights (scaled):'//
654 & 'WSCC= ',f10.6,' (SC-SC)'/
655 & 'WSCP= ',f10.6,' (SC-p)'/
656 & 'WELEC= ',f10.6,' (p-p electr)'/
657 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
658 & 'WBOND= ',f10.6,' (stretching)'/
659 & 'WANG= ',f10.6,' (bending)'/
660 & 'WSCLOC= ',f10.6,' (SC local)'/
661 & 'WTOR= ',f10.6,' (torsional)'/
662 & 'WTORD= ',f10.6,' (double torsional)'/
663 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
664 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
665 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
666 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
667 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
668 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
669 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
670 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
671 & 'WTURN6= ',f10.6,' (turns, 6th order)')
672 if(me.eq.king.or..not.out1file)
673 & write (iout,*) "Reference temperature for weights calculation:",
675 call reada(weightcard,"D0CM",d0cm,3.78d0)
676 call reada(weightcard,"AKCM",akcm,15.1d0)
677 call reada(weightcard,"AKTH",akth,11.0d0)
678 call reada(weightcard,"AKCT",akct,12.0d0)
679 call reada(weightcard,"V1SS",v1ss,-1.08d0)
680 call reada(weightcard,"V2SS",v2ss,7.61d0)
681 call reada(weightcard,"V3SS",v3ss,13.7d0)
682 call reada(weightcard,"EBR",ebr,-5.50D0)
683 call reada(weightcard,"ATRISS",atriss,0.301D0)
684 call reada(weightcard,"BTRISS",btriss,0.021D0)
685 call reada(weightcard,"CTRISS",ctriss,1.001D0)
686 call reada(weightcard,"DTRISS",dtriss,1.001D0)
687 write (iout,*) "ATRISS=", atriss
688 write (iout,*) "BTRISS=", btriss
689 write (iout,*) "CTRISS=", ctriss
690 write (iout,*) "DTRISS=", dtriss
691 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
693 dyn_ss_mask(i)=.false.
697 dyn_ssbond_ij(i,j)=1.0d300
700 call reada(weightcard,"HT",Ht,0.0D0)
702 ss_depth=ebr/wsc-0.25*eps(1,1)
703 Ht=Ht/wsc-0.25*eps(1,1)
704 akcm=akcm*wstrain/wsc
705 akth=akth*wstrain/wsc
706 akct=akct*wstrain/wsc
707 v1ss=v1ss*wstrain/wsc
708 v2ss=v2ss*wstrain/wsc
709 v3ss=v3ss*wstrain/wsc
711 if (wstrain.ne.0.0) then
712 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
718 if(me.eq.king.or..not.out1file) then
719 write (iout,*) "Parameters of the SS-bond potential:"
720 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
722 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
723 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
724 write (iout,*)" HT",Ht
725 print *,'indpdb=',indpdb,' pdbref=',pdbref
727 if (indpdb.gt.0 .or. pdbref) then
728 read(inp,'(a)') pdbfile
729 if(me.eq.king.or..not.out1file)
730 & write (iout,'(2a)') 'PDB data will be read from file ',
731 & pdbfile(:ilen(pdbfile))
732 open(ipdbin,file=pdbfile,status='old',err=33)
734 33 write (iout,'(a)') 'Error opening PDB file.'
737 c write (iout,*) 'Begin reading pdb data'
740 c write (iout,*) 'Finished reading pdb data'
742 if(me.eq.king.or..not.out1file)
743 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
744 & ' nstart_sup=',nstart_sup
746 itype_pdb(i)=itype(i)
750 nct=nstart_sup+nsup-1
751 call contact(.false.,ncont_ref,icont_ref,co)
754 if(me.eq.king.or..not.out1file)
755 & write(iout,*)'Adding sidechains'
759 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
762 do while (fail.and.nsi.le.maxsi)
763 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
766 if(fail) write(iout,*)'Adding sidechain failed for res ',
767 & i,' after ',nsi,' trials'
772 if (indpdb.eq.0) then
773 C Read sequence if not taken from the pdb file.
775 c print *,'nres=',nres
776 if (iscode.gt.0) then
777 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
779 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
781 C Convert sequence to numeric code
783 itype(i)=rescode(i,sequence(i),iscode)
785 C Assign initial virtual bond lengths
791 vbld(i+nres)=dsc(iabs(itype(i)))
792 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
793 c write (iout,*) "i",i," itype",itype(i),
794 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
798 c print '(20i4)',(itype(i),i=1,nres)
801 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
803 if (itype(i).eq.ntyp1) then
807 else if (iabs(itype(i+1)).ne.20) then
809 else if (iabs(itype(i)).ne.20) then
816 if(me.eq.king.or..not.out1file)then
817 write (iout,*) "ITEL"
819 write (iout,*) i,itype(i),itel(i)
821 print *,'Call Read_Bridge.'
824 C 8/13/98 Set limits to generating the dihedral angles
829 read (inp,*) ndih_constr
830 if (ndih_constr.gt.0) then
832 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
834 if(me.eq.king.or..not.out1file)then
836 & 'There are',ndih_constr,' constraints on phi angles.'
838 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
843 phi0(i)=deg2rad*phi0(i)
844 drange(i)=deg2rad*drange(i)
846 C if(me.eq.king.or..not.out1file)
847 C & write (iout,*) 'FTORS',ftors
850 phibound(1,ii) = phi0(i)-drange(i)
851 phibound(2,ii) = phi0(i)+drange(i)
858 write (iout,'(a)') 'Boundaries in phi angle sampling:'
860 write (iout,'(a3,i5,2f10.1)')
861 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
867 cd print *,'NNT=',NNT,' NCT=',NCT
868 if (itype(1).eq.ntyp1) nnt=2
869 if (itype(nres).eq.ntyp1) nct=nct-1
871 if(me.eq.king.or..not.out1file)
872 & write (iout,'(a,i3)') 'nsup=',nsup
874 if (nsup.le.(nct-nnt+1)) then
875 do i=0,nct-nnt+1-nsup
876 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
882 & 'Error - sequences to be superposed do not match.'
885 do i=0,nsup-(nct-nnt+1)
886 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
888 nstart_sup=nstart_sup+i
894 & 'Error - sequences to be superposed do not match.'
897 if (nsup.eq.0) nsup=nct-nnt
898 if (nstart_sup.eq.0) nstart_sup=nnt
899 if (nstart_seq.eq.0) nstart_seq=nnt
900 if(me.eq.king.or..not.out1file)
901 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
902 & ' nstart_seq=',nstart_seq
904 c--- Zscore rms -------
905 if (nz_start.eq.0) nz_start=nnt
906 if (nz_end.eq.0 .and. nsup.gt.0) then
908 else if (nz_end.eq.0) then
911 if(me.eq.king.or..not.out1file)then
912 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
913 write (iout,*) 'IZ_SC=',iz_sc
915 c----------------------
918 if (.not.pdbref) then
919 call read_angles(inp,*38)
921 38 write (iout,'(a)') 'Error reading reference structure.'
923 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
924 stop 'Error reading reference structure'
928 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
938 call contact(.true.,ncont_ref,icont_ref,co)
942 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
944 if (constr_dist.gt.0) call read_dist_constr
945 write (iout,*) "After read_dist_constr nhpb",nhpb
947 if(me.eq.king.or..not.out1file)
948 & write (iout,*) 'Contact order:',co
950 if(me.eq.king.or..not.out1file)
951 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
954 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
956 if(me.eq.king.or..not.out1file)
957 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
958 & icont_ref(1,i),' ',
959 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
963 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
964 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
965 & modecalc.ne.10) then
966 C If input structure hasn't been supplied from the PDB file read or generate
968 if (iranconf.eq.0 .and. .not. extconf) then
969 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
970 & write (iout,'(a)') 'Initial geometry will be read in.'
972 read(inp,'(8f10.5)',end=36,err=36)
973 & ((c(l,k),l=1,3),k=1,nres),
974 & ((c(l,k+nres),l=1,3),k=nnt,nct)
975 write (iout,*) "Exit READ_CART"
976 write (iout,'(8f10.5)')
977 & ((c(l,k),l=1,3),k=1,nres),
978 & ((c(l,k+nres),l=1,3),k=nnt,nct)
979 call int_from_cart1(.true.)
980 write (iout,*) "Finish INT_TO_CART"
983 dc(j,i)=c(j,i+1)-c(j,i)
984 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
988 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
990 dc(j,i+nres)=c(j,i+nres)-c(j,i)
991 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
997 call read_angles(inp,*36)
1000 36 write (iout,'(a)') 'Error reading angle file.'
1002 call mpi_finalize( MPI_COMM_WORLD,IERR )
1004 stop 'Error reading angle file.'
1006 else if (extconf) then
1007 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1008 & write (iout,'(a)') 'Extended chain initial geometry.'
1010 theta(i)=90d0*deg2rad
1013 phi(i)=180d0*deg2rad
1016 alph(i)=110d0*deg2rad
1019 omeg(i)=-120d0*deg2rad
1020 if (itype(i).le.0) omeg(i)=-omeg(i)
1023 if(me.eq.king.or..not.out1file)
1024 & write (iout,'(a)') 'Random-generated initial geometry.'
1028 if (me.eq.king .or. fg_rank.eq.0 .and. (
1029 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1033 call gen_rand_conf(itmp,*30)
1035 30 write (iout,*) 'Failed to generate random conformation',
1036 & ', itrial=',itrial
1037 write (*,*) 'Processor:',me,
1038 & ' Failed to generate random conformation',
1048 write (iout,'(a,i3,a)') 'Processor:',me,
1049 & ' error in generating random conformation.'
1050 write (*,'(a,i3,a)') 'Processor:',me,
1051 & ' error in generating random conformation.'
1054 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1060 call gen_rand_conf(itmp,*30)
1062 30 write (iout,*) 'Failed to generate random conformation',
1063 & ', itrial=',itrial
1064 write (*,*) 'Failed to generate random conformation',
1065 & ', itrial=',itrial
1067 write (iout,'(a,i3,a)') 'Processor:',me,
1068 & ' error in generating random conformation.'
1069 write (*,'(a,i3,a)') 'Processor:',me,
1070 & ' error in generating random conformation.'
1075 elseif (modecalc.eq.4) then
1076 read (inp,'(a)') intinname
1077 open (intin,file=intinname,status='old',err=333)
1078 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1079 & write (iout,'(a)') 'intinname',intinname
1080 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1082 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1084 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1086 stop 'Error opening angle file.'
1090 C Generate distance constraints, if the PDB structure is to be regularized.
1091 if (nthread.gt.0) then
1092 call read_threadbase
1095 if (me.eq.king .or. .not. out1file)
1097 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1098 write (iout,'(/a,i3,a)')
1099 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1100 write (iout,'(20i4)') (iss(i),i=1,ns)
1102 write(iout,*)"Running with dynamic disulfide-bond formation"
1104 write (iout,'(/a/)') 'Pre-formed links are:'
1110 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1111 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1117 if (ns.gt.0.and.dyn_ss) then
1121 forcon(i-nss)=forcon(i)
1128 dyn_ss_mask(iss(i))=.true.
1131 if (i2ndstr.gt.0) call secstrp2dihc
1132 c call geom_to_var(nvar,x)
1133 c call etotal(energia(0))
1134 c call enerprint(energia(0))
1135 c call briefout(0,etot)
1137 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1138 cd write (iout,'(a)') 'Variable list:'
1139 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1141 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1142 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1143 & 'Processor',myrank,': end reading molecular data.'
1148 c--------------------------------------------------------------------------
1149 logical function seq_comp(itypea,itypeb,length)
1151 integer length,itypea(length),itypeb(length)
1154 if (itypea(i).ne.itypeb(i)) then
1162 c-----------------------------------------------------------------------------
1163 subroutine read_bridge
1164 C Read information about disulfide bridges.
1165 implicit real*8 (a-h,o-z)
1166 include 'DIMENSIONS'
1170 include 'COMMON.IOUNITS'
1171 include 'COMMON.GEO'
1172 include 'COMMON.VAR'
1173 include 'COMMON.INTERACT'
1174 include 'COMMON.LOCAL'
1175 include 'COMMON.NAMES'
1176 include 'COMMON.CHAIN'
1177 include 'COMMON.FFIELD'
1178 include 'COMMON.SBRIDGE'
1179 include 'COMMON.HEADER'
1180 include 'COMMON.CONTROL'
1181 include 'COMMON.DBASE'
1182 include 'COMMON.THREAD'
1183 include 'COMMON.TIME1'
1184 include 'COMMON.SETUP'
1185 C Read bridging residues.
1186 read (inp,*) ns,(iss(i),i=1,ns)
1188 if(me.eq.king.or..not.out1file)
1189 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1190 C Check whether the specified bridging residues are cystines.
1192 if (itype(iss(i)).ne.1) then
1193 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1194 & 'Do you REALLY think that the residue ',
1195 & restyp(itype(iss(i))),i,
1196 & ' can form a disulfide bridge?!!!'
1197 write (*,'(2a,i3,a)')
1198 & 'Do you REALLY think that the residue ',
1199 & restyp(itype(iss(i))),i,
1200 & ' can form a disulfide bridge?!!!'
1202 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1207 C Read preformed bridges.
1209 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1211 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1214 C Check if the residues involved in bridges are in the specified list of
1215 C bridging residues.
1218 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1219 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1220 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1221 & ' contains residues present in other pairs.'
1222 write (*,'(a,i3,a)') 'Disulfide pair',i,
1223 & ' contains residues present in other pairs.'
1225 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1231 if (ihpb(i).eq.iss(j)) goto 10
1233 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1236 if (jhpb(i).eq.iss(j)) goto 20
1238 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1244 ihpb(i)=ihpb(i)+nres
1245 jhpb(i)=jhpb(i)+nres
1251 c----------------------------------------------------------------------------
1252 subroutine read_x(kanal,*)
1253 implicit real*8 (a-h,o-z)
1254 include 'DIMENSIONS'
1255 include 'COMMON.GEO'
1256 include 'COMMON.VAR'
1257 include 'COMMON.CHAIN'
1258 include 'COMMON.IOUNITS'
1259 include 'COMMON.CONTROL'
1260 include 'COMMON.LOCAL'
1261 include 'COMMON.INTERACT'
1262 c Read coordinates from input
1264 read(kanal,'(8f10.5)',end=10,err=10)
1265 & ((c(l,k),l=1,3),k=1,nres),
1266 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1269 c(j,2*nres)=c(j,nres)
1271 call int_from_cart1(.false.)
1274 dc(j,i)=c(j,i+1)-c(j,i)
1275 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1279 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1281 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1282 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1290 c----------------------------------------------------------------------------
1291 subroutine read_threadbase
1292 implicit real*8 (a-h,o-z)
1293 include 'DIMENSIONS'
1294 include 'COMMON.IOUNITS'
1295 include 'COMMON.GEO'
1296 include 'COMMON.VAR'
1297 include 'COMMON.INTERACT'
1298 include 'COMMON.LOCAL'
1299 include 'COMMON.NAMES'
1300 include 'COMMON.CHAIN'
1301 include 'COMMON.FFIELD'
1302 include 'COMMON.SBRIDGE'
1303 include 'COMMON.HEADER'
1304 include 'COMMON.CONTROL'
1305 include 'COMMON.DBASE'
1306 include 'COMMON.THREAD'
1307 include 'COMMON.TIME1'
1308 C Read pattern database for threading.
1309 read (icbase,*) nseq
1311 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1312 & nres_base(2,i),nres_base(3,i)
1313 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1315 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1316 c & nres_base(2,i),nres_base(3,i)
1317 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1321 if (weidis.eq.0.0D0) weidis=0.1D0
1330 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1331 write (iout,'(a,i5)') 'nexcl: ',nexcl
1332 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1335 c------------------------------------------------------------------------------
1336 subroutine setup_var
1337 implicit real*8 (a-h,o-z)
1338 include 'DIMENSIONS'
1339 include 'COMMON.IOUNITS'
1340 include 'COMMON.GEO'
1341 include 'COMMON.VAR'
1342 include 'COMMON.INTERACT'
1343 include 'COMMON.LOCAL'
1344 include 'COMMON.NAMES'
1345 include 'COMMON.CHAIN'
1346 include 'COMMON.FFIELD'
1347 include 'COMMON.SBRIDGE'
1348 include 'COMMON.HEADER'
1349 include 'COMMON.CONTROL'
1350 include 'COMMON.DBASE'
1351 include 'COMMON.THREAD'
1352 include 'COMMON.TIME1'
1353 C Set up variable list.
1359 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1361 ialph(i,1)=nvar+nside
1365 if (indphi.gt.0) then
1367 else if (indback.gt.0) then
1372 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1375 c----------------------------------------------------------------------------
1376 subroutine gen_dist_constr
1377 C Generate CA distance constraints.
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 dimension itype_pdb(maxres)
1395 common /pizda/ itype_pdb
1397 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1398 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1399 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1401 do i=nstart_sup,nstart_sup+nsup-1
1402 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1403 cd & ' seq_pdb', restyp(itype_pdb(i))
1404 do j=i+2,nstart_sup+nsup-1
1406 ihpb(nhpb)=i+nstart_seq-nstart_sup
1407 jhpb(nhpb)=j+nstart_seq-nstart_sup
1409 dhpb(nhpb)=dist(i,j)
1412 cd write (iout,'(a)') 'Distance constraints:'
1417 cd if (ii.gt.nres) then
1422 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1423 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1424 cd & dhpb(i),forcon(i)
1428 c----------------------------------------------------------------------------
1430 implicit real*8 (a-h,o-z)
1431 include 'DIMENSIONS'
1432 include 'COMMON.MAP'
1433 include 'COMMON.IOUNITS'
1434 character*3 angid(4) /'THE','PHI','ALP','OME'/
1435 character*80 mapcard,ucase
1437 read (inp,'(a)') mapcard
1438 mapcard=ucase(mapcard)
1439 if (index(mapcard,'PHI').gt.0) then
1441 else if (index(mapcard,'THE').gt.0) then
1443 else if (index(mapcard,'ALP').gt.0) then
1445 else if (index(mapcard,'OME').gt.0) then
1448 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1449 stop 'Error - illegal variable spec in MAP card.'
1451 call readi (mapcard,'RES1',res1(imap),0)
1452 call readi (mapcard,'RES2',res2(imap),0)
1453 if (res1(imap).eq.0) then
1454 res1(imap)=res2(imap)
1455 else if (res2(imap).eq.0) then
1456 res2(imap)=res1(imap)
1458 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1460 & 'Error - illegal definition of variable group in MAP.'
1461 stop 'Error - illegal definition of variable group in MAP.'
1463 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1464 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1465 call readi(mapcard,'NSTEP',nstep(imap),0)
1466 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1468 & 'Illegal boundary and/or step size specification in MAP.'
1469 stop 'Illegal boundary and/or step size specification in MAP.'
1474 c----------------------------------------------------------------------------
1476 implicit real*8 (a-h,o-z)
1477 include 'DIMENSIONS'
1478 include 'COMMON.IOUNITS'
1479 include 'COMMON.GEO'
1480 include 'COMMON.CSA'
1481 include 'COMMON.BANK'
1482 include 'COMMON.CONTROL'
1484 character*620 mcmcard
1485 call card_concat(mcmcard)
1487 call readi(mcmcard,'NCONF',nconf,50)
1488 call readi(mcmcard,'NADD',nadd,0)
1489 call readi(mcmcard,'JSTART',jstart,1)
1490 call readi(mcmcard,'JEND',jend,1)
1491 call readi(mcmcard,'NSTMAX',nstmax,500000)
1492 call readi(mcmcard,'N0',n0,1)
1493 call readi(mcmcard,'N1',n1,6)
1494 call readi(mcmcard,'N2',n2,4)
1495 call readi(mcmcard,'N3',n3,0)
1496 call readi(mcmcard,'N4',n4,0)
1497 call readi(mcmcard,'N5',n5,0)
1498 call readi(mcmcard,'N6',n6,10)
1499 call readi(mcmcard,'N7',n7,0)
1500 call readi(mcmcard,'N8',n8,0)
1501 call readi(mcmcard,'N9',n9,0)
1502 call readi(mcmcard,'N14',n14,0)
1503 call readi(mcmcard,'N15',n15,0)
1504 call readi(mcmcard,'N16',n16,0)
1505 call readi(mcmcard,'N17',n17,0)
1506 call readi(mcmcard,'N18',n18,0)
1508 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1510 call readi(mcmcard,'NDIFF',ndiff,2)
1511 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1512 call readi(mcmcard,'IS1',is1,1)
1513 call readi(mcmcard,'IS2',is2,8)
1514 call readi(mcmcard,'NRAN0',nran0,4)
1515 call readi(mcmcard,'NRAN1',nran1,2)
1516 call readi(mcmcard,'IRR',irr,1)
1517 call readi(mcmcard,'NSEED',nseed,20)
1518 call readi(mcmcard,'NTOTAL',ntotal,10000)
1519 call reada(mcmcard,'CUT1',cut1,2.0d0)
1520 call reada(mcmcard,'CUT2',cut2,5.0d0)
1521 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1522 call readi(mcmcard,'ICMAX',icmax,3)
1523 call readi(mcmcard,'IRESTART',irestart,0)
1524 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1527 call reada(mcmcard,'DELE',dele,20.0d0)
1528 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1529 call readi(mcmcard,'IREF',iref,0)
1530 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1531 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1532 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1533 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1534 write (iout,*) "NCONF_IN",nconf_in
1537 c----------------------------------------------------------------------------
1538 cfmc subroutine mcmfread
1539 cfmc implicit real*8 (a-h,o-z)
1540 cfmc include 'DIMENSIONS'
1541 cfmc include 'COMMON.MCMF'
1542 cfmc include 'COMMON.IOUNITS'
1543 cfmc include 'COMMON.GEO'
1544 cfmc character*80 ucase
1545 cfmc character*620 mcmcard
1546 cfmc call card_concat(mcmcard)
1548 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1549 cfmc write(iout,*)'MAXRANT=',maxrant
1550 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1551 cfmc write(iout,*)'MAXFAM=',maxfam
1552 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1553 cfmc write(iout,*)'NNET1=',nnet1
1554 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1555 cfmc write(iout,*)'NNET2=',nnet2
1556 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1557 cfmc write(iout,*)'NNET3=',nnet3
1558 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1559 cfmc write(iout,*)'ILASTT=',ilastt
1560 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1561 cfmc write(iout,*)'MAXSTR=',maxstr
1562 cfmc maxstr_f=maxstr/maxfam
1563 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1564 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1565 cfmc write(iout,*)'NMCMF=',nmcmf
1566 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1567 cfmc write(iout,*)'IFOCUS=',ifocus
1568 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1569 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1570 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1571 cfmc write(iout,*)'INTPRT=',intprt
1572 cfmc call readi(mcmcard,'IPRT',iprt,100)
1573 cfmc write(iout,*)'IPRT=',iprt
1574 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1575 cfmc write(iout,*)'IMAXTR=',imaxtr
1576 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1577 cfmc write(iout,*)'MAXEVEN=',maxeven
1578 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1579 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1580 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1581 cfmc write(iout,*)'INIMIN=',inimin
1582 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1583 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1584 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1585 cfmc write(iout,*)'NTHREAD=',nthread
1586 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1587 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1588 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1589 cfmc write(iout,*)'MAXPERT=',maxpert
1590 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1591 cfmc write(iout,*)'IRMSD=',irmsd
1592 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1593 cfmc write(iout,*)'DENEMIN=',denemin
1594 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1595 cfmc write(iout,*)'RCUT1S=',rcut1s
1596 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1597 cfmc write(iout,*)'RCUT1E=',rcut1e
1598 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1599 cfmc write(iout,*)'RCUT2S=',rcut2s
1600 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1601 cfmc write(iout,*)'RCUT2E=',rcut2e
1602 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1603 cfmc write(iout,*)'DPERT1=',d_pert1
1604 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1605 cfmc write(iout,*)'DPERT1A=',d_pert1a
1606 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1607 cfmc write(iout,*)'DPERT2=',d_pert2
1608 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1609 cfmc write(iout,*)'DPERT2A=',d_pert2a
1610 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1611 cfmc write(iout,*)'DPERT2B=',d_pert2b
1612 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1613 cfmc write(iout,*)'DPERT2C=',d_pert2c
1614 cfmc d_pert1=deg2rad*d_pert1
1615 cfmc d_pert1a=deg2rad*d_pert1a
1616 cfmc d_pert2=deg2rad*d_pert2
1617 cfmc d_pert2a=deg2rad*d_pert2a
1618 cfmc d_pert2b=deg2rad*d_pert2b
1619 cfmc d_pert2c=deg2rad*d_pert2c
1620 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1621 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1622 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1623 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1624 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1625 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1626 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1627 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1628 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1629 cfmc write(iout,*)'RCUTINI=',rcutini
1630 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1631 cfmc write(iout,*)'GRAT=',grat
1632 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1633 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1637 c----------------------------------------------------------------------------
1639 implicit real*8 (a-h,o-z)
1640 include 'DIMENSIONS'
1641 include 'COMMON.MCM'
1642 include 'COMMON.MCE'
1643 include 'COMMON.IOUNITS'
1645 character*320 mcmcard
1646 call card_concat(mcmcard)
1647 call readi(mcmcard,'MAXACC',maxacc,100)
1648 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1649 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1650 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1651 call readi(mcmcard,'MAXREPM',maxrepm,200)
1652 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1653 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1654 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1655 call reada(mcmcard,'E_UP',e_up,5.0D0)
1656 call reada(mcmcard,'DELTE',delte,0.1D0)
1657 call readi(mcmcard,'NSWEEP',nsweep,5)
1658 call readi(mcmcard,'NSTEPH',nsteph,0)
1659 call readi(mcmcard,'NSTEPC',nstepc,0)
1660 call reada(mcmcard,'TMIN',tmin,298.0D0)
1661 call reada(mcmcard,'TMAX',tmax,298.0D0)
1662 call readi(mcmcard,'NWINDOW',nwindow,0)
1663 call readi(mcmcard,'PRINT_MC',print_mc,0)
1664 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1665 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1666 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1667 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1668 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1669 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1670 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1671 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1672 if (nwindow.gt.0) then
1673 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1675 winlen(i)=winend(i)-winstart(i)+1
1678 if (tmax.lt.tmin) tmax=tmin
1679 if (tmax.eq.tmin) then
1683 if (nstepc.gt.0 .and. nsteph.gt.0) then
1684 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1685 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1687 C Probabilities of different move types
1688 sumpro_type(0)=0.0D0
1689 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1690 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1691 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1692 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1693 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1694 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1695 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1697 print *,'i',i,' sumprotype',sumpro_type(i)
1698 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1699 print *,'i',i,' sumprotype',sumpro_type(i)
1703 c----------------------------------------------------------------------------
1704 subroutine read_minim
1705 implicit real*8 (a-h,o-z)
1706 include 'DIMENSIONS'
1707 include 'COMMON.MINIM'
1708 include 'COMMON.IOUNITS'
1710 character*320 minimcard
1711 call card_concat(minimcard)
1712 call readi(minimcard,'MAXMIN',maxmin,2000)
1713 call readi(minimcard,'MAXFUN',maxfun,5000)
1714 call readi(minimcard,'MINMIN',minmin,maxmin)
1715 call readi(minimcard,'MINFUN',minfun,maxmin)
1716 call reada(minimcard,'TOLF',tolf,1.0D-2)
1717 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1718 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1719 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1720 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1721 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1722 & 'Options in energy minimization:'
1723 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1724 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1725 & 'MinMin:',MinMin,' MinFun:',MinFun,
1726 & ' TolF:',TolF,' RTolF:',RTolF
1729 c----------------------------------------------------------------------------
1730 subroutine read_angles(kanal,*)
1731 implicit real*8 (a-h,o-z)
1732 include 'DIMENSIONS'
1733 include 'COMMON.GEO'
1734 include 'COMMON.VAR'
1735 include 'COMMON.CHAIN'
1736 include 'COMMON.IOUNITS'
1737 include 'COMMON.CONTROL'
1738 c Read angles from input
1740 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1741 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1742 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1743 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1746 c 9/7/01 avoid 180 deg valence angle
1747 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1749 theta(i)=deg2rad*theta(i)
1750 phi(i)=deg2rad*phi(i)
1751 alph(i)=deg2rad*alph(i)
1752 omeg(i)=deg2rad*omeg(i)
1757 c----------------------------------------------------------------------------
1758 subroutine reada(rekord,lancuch,wartosc,default)
1760 character*(*) rekord,lancuch
1761 double precision wartosc,default
1764 iread=index(rekord,lancuch)
1765 if (iread.eq.0) then
1769 iread=iread+ilen(lancuch)+1
1770 read (rekord(iread:),*,err=10,end=10) wartosc
1775 c----------------------------------------------------------------------------
1776 subroutine readi(rekord,lancuch,wartosc,default)
1778 character*(*) rekord,lancuch
1779 integer wartosc,default
1782 iread=index(rekord,lancuch)
1783 if (iread.eq.0) then
1787 iread=iread+ilen(lancuch)+1
1788 read (rekord(iread:),*,err=10,end=10) wartosc
1793 c----------------------------------------------------------------------------
1794 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1797 integer tablica(dim),default
1798 character*(*) rekord,lancuch
1805 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1806 if (iread.eq.0) return
1807 iread=iread+ilen(lancuch)+1
1808 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1811 c----------------------------------------------------------------------------
1812 subroutine multreada(rekord,lancuch,tablica,dim,default)
1815 double precision tablica(dim),default
1816 character*(*) rekord,lancuch
1823 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1824 if (iread.eq.0) return
1825 iread=iread+ilen(lancuch)+1
1826 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1829 c----------------------------------------------------------------------------
1830 subroutine openunits
1831 implicit real*8 (a-h,o-z)
1832 include 'DIMENSIONS'
1835 character*16 form,nodename
1838 include 'COMMON.SETUP'
1839 include 'COMMON.IOUNITS'
1841 include 'COMMON.CONTROL'
1842 integer lenpre,lenpot,ilen,lentmp
1844 character*3 out1file_text,ucase
1847 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1848 call getenv_loc("PREFIX",prefix)
1850 call getenv_loc("POT",pot)
1851 call getenv_loc("DIRTMP",tmpdir)
1852 call getenv_loc("CURDIR",curdir)
1853 call getenv_loc("OUT1FILE",out1file_text)
1854 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1855 out1file_text=ucase(out1file_text)
1856 if (out1file_text(1:1).eq."Y") then
1859 out1file=fg_rank.gt.0
1864 if (lentmp.gt.0) then
1865 write (*,'(80(1h!))')
1866 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1867 write (*,'(80(1h!))')
1868 write (*,*)"All output files will be on node /tmp directory."
1870 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1871 if (me.eq.king) then
1872 write (*,*) "The master node is ",nodename
1873 else if (fg_rank.eq.0) then
1874 write (*,*) "I am the CG slave node ",nodename
1876 write (*,*) "I am the FG slave node ",nodename
1879 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1880 lenpre = lentmp+lenpre+1
1882 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1883 C Get the names and open the input files
1884 #if defined(WINIFL) || defined(WINPGI)
1885 open(1,file=pref_orig(:ilen(pref_orig))//
1886 & '.inp',status='old',readonly,shared)
1887 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1888 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1889 C Get parameter filenames and open the parameter files.
1890 call getenv_loc('BONDPAR',bondname)
1891 open (ibond,file=bondname,status='old',readonly,shared)
1892 call getenv_loc('THETPAR',thetname)
1893 open (ithep,file=thetname,status='old',readonly,shared)
1894 call getenv_loc('ROTPAR',rotname)
1895 open (irotam,file=rotname,status='old',readonly,shared)
1896 call getenv_loc('TORPAR',torname)
1897 open (itorp,file=torname,status='old',readonly,shared)
1898 call getenv_loc('TORDPAR',tordname)
1899 open (itordp,file=tordname,status='old',readonly,shared)
1900 call getenv_loc('FOURIER',fouriername)
1901 open (ifourier,file=fouriername,status='old',readonly,shared)
1902 call getenv_loc('ELEPAR',elename)
1903 open (ielep,file=elename,status='old',readonly,shared)
1904 call getenv_loc('SIDEPAR',sidename)
1905 open (isidep,file=sidename,status='old',readonly,shared)
1906 #elif (defined CRAY) || (defined AIX)
1907 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1909 c print *,"Processor",myrank," opened file 1"
1910 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1911 c print *,"Processor",myrank," opened file 9"
1912 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1913 C Get parameter filenames and open the parameter files.
1914 call getenv_loc('BONDPAR',bondname)
1915 open (ibond,file=bondname,status='old',action='read')
1916 c print *,"Processor",myrank," opened file IBOND"
1917 call getenv_loc('THETPAR',thetname)
1918 open (ithep,file=thetname,status='old',action='read')
1919 c print *,"Processor",myrank," opened file ITHEP"
1920 call getenv_loc('ROTPAR',rotname)
1921 open (irotam,file=rotname,status='old',action='read')
1922 c print *,"Processor",myrank," opened file IROTAM"
1923 call getenv_loc('TORPAR',torname)
1924 open (itorp,file=torname,status='old',action='read')
1925 c print *,"Processor",myrank," opened file ITORP"
1926 call getenv_loc('TORDPAR',tordname)
1927 open (itordp,file=tordname,status='old',action='read')
1928 c print *,"Processor",myrank," opened file ITORDP"
1929 call getenv_loc('SCCORPAR',sccorname)
1930 open (isccor,file=sccorname,status='old',action='read')
1931 c print *,"Processor",myrank," opened file ISCCOR"
1932 call getenv_loc('FOURIER',fouriername)
1933 open (ifourier,file=fouriername,status='old',action='read')
1934 c print *,"Processor",myrank," opened file IFOURIER"
1935 call getenv_loc('ELEPAR',elename)
1936 open (ielep,file=elename,status='old',action='read')
1937 c print *,"Processor",myrank," opened file IELEP"
1938 call getenv_loc('SIDEPAR',sidename)
1939 open (isidep,file=sidename,status='old',action='read')
1940 c print *,"Processor",myrank," opened file ISIDEP"
1941 c print *,"Processor",myrank," opened parameter files"
1943 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1944 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1945 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1946 C Get parameter filenames and open the parameter files.
1947 call getenv_loc('BONDPAR',bondname)
1948 open (ibond,file=bondname,status='old')
1949 call getenv_loc('THETPAR',thetname)
1950 open (ithep,file=thetname,status='old')
1951 call getenv_loc('ROTPAR',rotname)
1952 open (irotam,file=rotname,status='old')
1953 call getenv_loc('TORPAR',torname)
1954 open (itorp,file=torname,status='old')
1955 call getenv_loc('TORDPAR',tordname)
1956 open (itordp,file=tordname,status='old')
1957 call getenv_loc('SCCORPAR',sccorname)
1958 open (isccor,file=sccorname,status='old')
1959 call getenv_loc('FOURIER',fouriername)
1960 open (ifourier,file=fouriername,status='old')
1961 call getenv_loc('ELEPAR',elename)
1962 open (ielep,file=elename,status='old')
1963 call getenv_loc('SIDEPAR',sidename)
1964 open (isidep,file=sidename,status='old')
1966 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1968 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1969 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1970 C Get parameter filenames and open the parameter files.
1971 call getenv_loc('BONDPAR',bondname)
1972 open (ibond,file=bondname,status='old',readonly)
1973 call getenv_loc('THETPAR',thetname)
1974 open (ithep,file=thetname,status='old',readonly)
1975 call getenv_loc('ROTPAR',rotname)
1976 open (irotam,file=rotname,status='old',readonly)
1977 call getenv_loc('TORPAR',torname)
1978 open (itorp,file=torname,status='old',readonly)
1979 call getenv_loc('TORDPAR',tordname)
1980 open (itordp,file=tordname,status='old',readonly)
1981 call getenv_loc('SCCORPAR',sccorname)
1982 open (isccor,file=sccorname,status='old',readonly)
1984 call getenv_loc('THETPARPDB',thetname_pdb)
1985 print *,"thetname_pdb ",thetname_pdb
1986 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1987 print *,ithep_pdb," opened"
1989 call getenv_loc('FOURIER',fouriername)
1990 open (ifourier,file=fouriername,status='old',readonly)
1991 call getenv_loc('ELEPAR',elename)
1992 open (ielep,file=elename,status='old',readonly)
1993 call getenv_loc('SIDEPAR',sidename)
1994 open (isidep,file=sidename,status='old',readonly)
1996 call getenv_loc('ROTPARPDB',rotname_pdb)
1997 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2002 C 8/9/01 In the newest version SCp interaction constants are read from a file
2003 C Use -DOLDSCP to use hard-coded constants instead.
2005 call getenv_loc('SCPPAR',scpname)
2006 #if defined(WINIFL) || defined(WINPGI)
2007 open (iscpp,file=scpname,status='old',readonly,shared)
2008 #elif (defined CRAY) || (defined AIX)
2009 open (iscpp,file=scpname,status='old',action='read')
2011 open (iscpp,file=scpname,status='old')
2013 open (iscpp,file=scpname,status='old',readonly)
2016 call getenv_loc('PATTERN',patname)
2017 #if defined(WINIFL) || defined(WINPGI)
2018 open (icbase,file=patname,status='old',readonly,shared)
2019 #elif (defined CRAY) || (defined AIX)
2020 open (icbase,file=patname,status='old',action='read')
2022 open (icbase,file=patname,status='old')
2024 open (icbase,file=patname,status='old',readonly)
2027 C Open output file only for CG processes
2028 c print *,"Processor",myrank," fg_rank",fg_rank
2029 if (fg_rank.eq.0) then
2031 if (nodes.eq.1) then
2034 npos = dlog10(dfloat(nodes-1))+1
2036 if (npos.lt.3) npos=3
2037 write (liczba,'(i1)') npos
2038 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2040 write (liczba,form) me
2041 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2042 & liczba(:ilen(liczba))
2043 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2045 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2047 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2048 & liczba(:ilen(liczba))//'.mol2'
2049 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2050 & liczba(:ilen(liczba))//'.stat'
2052 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2053 & //liczba(:ilen(liczba))//'.stat')
2054 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2057 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2058 & liczba(:ilen(liczba))//'.const'
2063 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2064 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2065 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2066 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2067 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2069 & call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
2071 rest2name=prefix(:ilen(prefix))//'.rst'
2073 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2076 #if defined(AIX) || defined(PGI)
2077 if (me.eq.king .or. .not. out1file)
2078 & open(iout,file=outname,status='unknown')
2080 if (fg_rank.gt.0) then
2081 write (liczba,'(i3.3)') myrank/nfgtasks
2082 write (ll,'(bz,i3.3)') fg_rank
2083 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2088 open(igeom,file=intname,status='unknown',position='append')
2089 open(ipdb,file=pdbname,status='unknown')
2090 open(imol2,file=mol2name,status='unknown')
2091 open(istat,file=statname,status='unknown',position='append')
2093 c1out open(iout,file=outname,status='unknown')
2096 if (me.eq.king .or. .not.out1file)
2097 & open(iout,file=outname,status='unknown')
2099 if (fg_rank.gt.0) then
2100 write (liczba,'(i3.3)') myrank/nfgtasks
2101 write (ll,'(bz,i3.3)') fg_rank
2102 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2107 open(igeom,file=intname,status='unknown',access='append')
2108 open(ipdb,file=pdbname,status='unknown')
2109 open(imol2,file=mol2name,status='unknown')
2110 open(istat,file=statname,status='unknown',access='append')
2112 c1out open(iout,file=outname,status='unknown')
2115 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2116 csa_seed=prefix(:lenpre)//'.CSA.seed'
2117 csa_history=prefix(:lenpre)//'.CSA.history'
2118 csa_bank=prefix(:lenpre)//'.CSA.bank'
2119 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2120 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2121 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2122 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2123 csa_int=prefix(:lenpre)//'.int'
2124 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2125 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2126 csa_in=prefix(:lenpre)//'.CSA.in'
2127 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2130 write (iout,'(80(1h-))')
2131 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2132 write (iout,'(80(1h-))')
2133 write (iout,*) "Input file : ",
2134 & pref_orig(:ilen(pref_orig))//'.inp'
2135 write (iout,*) "Output file : ",
2136 & outname(:ilen(outname))
2138 write (iout,*) "Sidechain potential file : ",
2139 & sidename(:ilen(sidename))
2141 write (iout,*) "SCp potential file : ",
2142 & scpname(:ilen(scpname))
2144 write (iout,*) "Electrostatic potential file : ",
2145 & elename(:ilen(elename))
2146 write (iout,*) "Cumulant coefficient file : ",
2147 & fouriername(:ilen(fouriername))
2148 write (iout,*) "Torsional parameter file : ",
2149 & torname(:ilen(torname))
2150 write (iout,*) "Double torsional parameter file : ",
2151 & tordname(:ilen(tordname))
2152 write (iout,*) "SCCOR parameter file : ",
2153 & sccorname(:ilen(sccorname))
2154 write (iout,*) "Bond & inertia constant file : ",
2155 & bondname(:ilen(bondname))
2156 write (iout,*) "Bending parameter file : ",
2157 & thetname(:ilen(thetname))
2158 write (iout,*) "Rotamer parameter file : ",
2159 & rotname(:ilen(rotname))
2160 write (iout,*) "Thetpdb parameter file : ",
2161 & thetname_pdb(:ilen(thetname_pdb))
2162 write (iout,*) "Threading database : ",
2163 & patname(:ilen(patname))
2165 &write (iout,*)" DIRTMP : ",
2167 write (iout,'(80(1h-))')
2171 c----------------------------------------------------------------------------
2172 subroutine card_concat(card)
2173 implicit real*8 (a-h,o-z)
2174 include 'DIMENSIONS'
2175 include 'COMMON.IOUNITS'
2177 character*80 karta,ucase
2179 read (inp,'(a)') karta
2182 do while (karta(80:80).eq.'&')
2183 card=card(:ilen(card)+1)//karta(:79)
2184 read (inp,'(a)') karta
2187 card=card(:ilen(card)+1)//karta
2190 c----------------------------------------------------------------------------------
2192 implicit real*8 (a-h,o-z)
2193 include 'DIMENSIONS'
2194 include 'COMMON.CHAIN'
2195 include 'COMMON.IOUNITS'
2197 open(irest2,file=rest2name,status='unknown')
2198 read(irest2,*) totT,EK,potE,totE,t_bath
2200 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2203 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2206 read (irest2,*) iset
2211 c---------------------------------------------------------------------------------
2212 subroutine read_fragments
2213 implicit real*8 (a-h,o-z)
2214 include 'DIMENSIONS'
2218 include 'COMMON.SETUP'
2219 include 'COMMON.CHAIN'
2220 include 'COMMON.IOUNITS'
2222 include 'COMMON.CONTROL'
2223 read(inp,*) nset,nfrag,npair,nfrag_back
2224 if(me.eq.king.or..not.out1file)
2225 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2226 & " nfrag_back",nfrag_back
2228 read(inp,*) mset(iset)
2230 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2232 if(me.eq.king.or..not.out1file)
2233 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2234 & ifrag(2,i,iset), qinfrag(i,iset)
2237 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2239 if(me.eq.king.or..not.out1file)
2240 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2241 & ipair(2,i,iset), qinpair(i,iset)
2244 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2245 & wfrag_back(3,i,iset),
2246 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2247 if(me.eq.king.or..not.out1file)
2248 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2249 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2254 c-------------------------------------------------------------------------------
2255 subroutine read_dist_constr
2256 implicit real*8 (a-h,o-z)
2257 include 'DIMENSIONS'
2261 include 'COMMON.SETUP'
2262 include 'COMMON.CONTROL'
2263 include 'COMMON.CHAIN'
2264 include 'COMMON.IOUNITS'
2265 include 'COMMON.SBRIDGE'
2266 integer ifrag_(2,100),ipair_(2,100)
2267 double precision wfrag_(100),wpair_(100)
2268 character*500 controlcard
2270 write (iout,*) "Calling read_dist_constr"
2271 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2273 call card_concat(controlcard)
2274 call readi(controlcard,"NFRAG",nfrag_,0)
2275 call readi(controlcard,"NPAIR",npair_,0)
2276 call readi(controlcard,"NDIST",ndist_,0)
2277 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2278 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2279 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2280 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2281 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2282 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2283 c write (iout,*) "IFRAG"
2285 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2287 c write (iout,*) "IPAIR"
2289 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2293 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2294 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2295 & ifrag_(2,i)=nstart_sup+nsup-1
2296 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2298 if (wfrag_(i).gt.0.0d0) then
2299 do j=ifrag_(1,i),ifrag_(2,i)-1
2300 do k=j+1,ifrag_(2,i)
2301 c write (iout,*) "j",j," k",k
2303 if (constr_dist.eq.1) then
2308 forcon(nhpb)=wfrag_(i)
2309 else if (constr_dist.eq.2) then
2310 if (ddjk.le.dist_cut) then
2315 forcon(nhpb)=wfrag_(i)
2322 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2325 if (.not.out1file .or. me.eq.king)
2326 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2327 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2329 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2330 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2337 if (wpair_(i).gt.0.0d0) then
2345 do j=ifrag_(1,ii),ifrag_(2,ii)
2346 do k=ifrag_(1,jj),ifrag_(2,jj)
2350 forcon(nhpb)=wpair_(i)
2351 dhpb(nhpb)=dist(j,k)
2353 if (.not.out1file .or. me.eq.king)
2354 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2355 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2357 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2358 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2366 if (constr_dist.eq.11) then
2367 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2368 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2369 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2372 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2373 & ibecarb(i),forcon(nhpb+1)
2375 if (forcon(nhpb+1).gt.0.0d0) then
2377 if (ibecarb(i).gt.0) then
2378 ihpb(i)=ihpb(i)+nres
2379 jhpb(i)=jhpb(i)+nres
2381 if (dhpb(nhpb).eq.0.0d0)
2382 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2384 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2385 C if (forcon(nhpb+1).gt.0.0d0) then
2387 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2389 if (.not.out1file .or. me.eq.king)
2390 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2391 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2393 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2394 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2401 c-------------------------------------------------------------------------------
2403 subroutine flush(iu)
2408 subroutine flush(iu)
2413 c------------------------------------------------------------------------------
2414 subroutine copy_to_tmp(source)
2415 include "DIMENSIONS"
2416 include "COMMON.IOUNITS"
2417 character*(*) source
2418 character* 256 tmpfile
2422 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2423 inquire(file=tmpfile,exist=ex)
2425 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2426 & " to temporary directory..."
2427 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2428 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2432 c------------------------------------------------------------------------------
2433 subroutine move_from_tmp(source)
2434 include "DIMENSIONS"
2435 include "COMMON.IOUNITS"
2436 character*(*) source
2439 write (*,*) "Moving ",source(:ilen(source)),
2440 & " from temporary directory to working directory"
2441 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2442 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2445 c------------------------------------------------------------------------------
2446 subroutine random_init(seed)
2448 C Initialize random number generator
2450 implicit real*8 (a-h,o-z)
2451 include 'DIMENSIONS'
2454 logical OKRandom, prng_restart
2456 integer iseed_array(4)
2458 include 'COMMON.IOUNITS'
2459 include 'COMMON.TIME1'
2460 include 'COMMON.THREAD'
2461 include 'COMMON.SBRIDGE'
2462 include 'COMMON.CONTROL'
2463 include 'COMMON.MCM'
2464 include 'COMMON.MAP'
2465 include 'COMMON.HEADER'
2466 include 'COMMON.CSA'
2467 include 'COMMON.CHAIN'
2468 include 'COMMON.MUCA'
2470 include 'COMMON.FFIELD'
2471 include 'COMMON.SETUP'
2472 iseed=-dint(dabs(seed))
2473 if (iseed.eq.0) then
2474 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2475 & 'Random seed undefined. The program will stop.'
2476 write (*,'(/80(1h*)/20x,a/80(1h*))')
2477 & 'Random seed undefined. The program will stop.'
2479 call mpi_finalize(mpi_comm_world,ierr)
2481 stop 'Bad random seed.'
2484 if (fg_rank.eq.0) then
2488 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2489 OKRandom = prng_restart(me,iseed)
2492 tmp=65536.0d0**(4-i)
2493 iseed_array(i) = dint(seed/tmp)
2494 seed=seed-iseed_array(i)*tmp
2497 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2498 & (iseed_array(i),i=1,4)
2499 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2500 & (iseed_array(i),i=1,4)
2501 OKRandom = prng_restart(me,iseed_array)
2504 c r1 = prng_next(me)
2505 r1=ran_number(0.0D0,1.0D0)
2507 & write (iout,*) 'ran_num',r1
2508 if (r1.lt.0.0d0) OKRandom=.false.
2510 if (.not.OKRandom) then
2511 write (iout,*) 'PRNG IS NOT WORKING!!!'
2512 print *,'PRNG IS NOT WORKING!!!'
2515 call mpi_abort(mpi_comm_world,error_msg,ierr)
2518 write (iout,*) 'too many processors for parallel prng'
2519 write (*,*) 'too many processors for parallel prng'
2527 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)