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 C DISTCHAINMAX become obsolete for periodic boundry condition
217 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
218 C Reading the dimensions of box in x,y,z coordinates
219 call reada(controlcard,'BOXX',boxxsize,100.0d0)
220 call reada(controlcard,'BOXY',boxysize,100.0d0)
221 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
223 if (me.eq.king .or. .not.out1file )
224 & write (iout,*) "DISTCHAINMAX",distchainmax
226 if(me.eq.king.or..not.out1file)
227 & write (iout,'(2a)') diagmeth(kdiag),
228 & ' routine used to diagonalize matrices.'
231 c--------------------------------------------------------------------------
232 subroutine read_REMDpar
236 implicit real*8 (a-h,o-z)
238 include 'COMMON.IOUNITS'
239 include 'COMMON.TIME1'
242 include 'COMMON.LANGEVIN'
244 include 'COMMON.LANGEVIN.lang0'
246 include 'COMMON.INTERACT'
247 include 'COMMON.NAMES'
249 include 'COMMON.REMD'
250 include 'COMMON.CONTROL'
251 include 'COMMON.SETUP'
253 character*320 controlcard
254 character*3200 controlcard1
255 integer iremd_m_total
257 if(me.eq.king.or..not.out1file)
258 & write (iout,*) "REMD setup"
260 call card_concat(controlcard)
261 call readi(controlcard,"NREP",nrep,3)
262 call readi(controlcard,"NSTEX",nstex,1000)
263 call reada(controlcard,"RETMIN",retmin,10.0d0)
264 call reada(controlcard,"RETMAX",retmax,1000.0d0)
265 mremdsync=(index(controlcard,'SYNC').gt.0)
266 call readi(controlcard,"NSYN",i_sync_step,100)
267 restart1file=(index(controlcard,'REST1FILE').gt.0)
268 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
269 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
270 if(max_cache_traj_use.gt.max_cache_traj)
271 & max_cache_traj_use=max_cache_traj
272 if(me.eq.king.or..not.out1file) then
273 cd if (traj1file) then
274 crc caching is in testing - NTWX is not ignored
275 cd write (iout,*) "NTWX value is ignored"
276 cd write (iout,*) " trajectory is stored to one file by master"
277 cd write (iout,*) " before exchange at NSTEX intervals"
279 write (iout,*) "NREP= ",nrep
280 write (iout,*) "NSTEX= ",nstex
281 write (iout,*) "SYNC= ",mremdsync
282 write (iout,*) "NSYN= ",i_sync_step
283 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
286 if (index(controlcard,'TLIST').gt.0) then
288 call card_concat(controlcard1)
289 read(controlcard1,*) (remd_t(i),i=1,nrep)
290 if(me.eq.king.or..not.out1file)
291 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
294 if (index(controlcard,'MLIST').gt.0) then
296 call card_concat(controlcard1)
297 read(controlcard1,*) (remd_m(i),i=1,nrep)
298 if(me.eq.king.or..not.out1file) then
299 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
302 iremd_m_total=iremd_m_total+remd_m(i)
304 write (iout,*) 'Total number of replicas ',iremd_m_total
307 if(me.eq.king.or..not.out1file)
308 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
311 c--------------------------------------------------------------------------
312 subroutine read_MDpar
316 implicit real*8 (a-h,o-z)
318 include 'COMMON.IOUNITS'
319 include 'COMMON.TIME1'
322 include 'COMMON.LANGEVIN'
324 include 'COMMON.LANGEVIN.lang0'
326 include 'COMMON.INTERACT'
327 include 'COMMON.NAMES'
329 include 'COMMON.SETUP'
330 include 'COMMON.CONTROL'
331 include 'COMMON.SPLITELE'
333 character*320 controlcard
335 call card_concat(controlcard)
336 call readi(controlcard,"NSTEP",n_timestep,1000000)
337 call readi(controlcard,"NTWE",ntwe,100)
338 call readi(controlcard,"NTWX",ntwx,1000)
339 call reada(controlcard,"DT",d_time,1.0d-1)
340 call reada(controlcard,"DVMAX",dvmax,2.0d1)
341 call reada(controlcard,"DAMAX",damax,1.0d1)
342 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
343 call readi(controlcard,"LANG",lang,0)
344 RESPA = index(controlcard,"RESPA") .gt. 0
345 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
346 ntime_split0=ntime_split
347 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
348 ntime_split0=ntime_split
349 call reada(controlcard,"R_CUT",r_cut,2.0d0)
350 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
351 rest = index(controlcard,"REST").gt.0
352 tbf = index(controlcard,"TBF").gt.0
353 usampl = index(controlcard,"USAMPL").gt.0
354 mdpdb = index(controlcard,"MDPDB").gt.0
355 call reada(controlcard,"T_BATH",t_bath,300.0d0)
356 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
357 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
358 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
359 if (count_reset_moment.eq.0) count_reset_moment=1000000000
360 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
361 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
362 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
363 if (count_reset_vel.eq.0) count_reset_vel=1000000000
364 large = index(controlcard,"LARGE").gt.0
365 print_compon = index(controlcard,"PRINT_COMPON").gt.0
366 rattle = index(controlcard,"RATTLE").gt.0
367 c if performing umbrella sampling, fragments constrained are read from the fragment file
373 if(me.eq.king.or..not.out1file) then
375 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
377 write (iout,'(a)') "The units are:"
378 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
379 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
380 & " acceleration: angstrom/(48.9 fs)**2"
381 write (iout,'(a)') "energy: kcal/mol, temperature: K"
383 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
384 write (iout,'(a60,f10.5,a)')
385 & "Initial time step of numerical integration:",d_time,
387 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
389 write (iout,'(2a,i4,a)')
390 & "A-MTS algorithm used; initial time step for fast-varying",
391 & " short-range forces split into",ntime_split," steps."
392 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
393 & r_cut," lambda",rlamb
395 write (iout,'(2a,f10.5)')
396 & "Maximum acceleration threshold to reduce the time step",
397 & "/increase split number:",damax
398 write (iout,'(2a,f10.5)')
399 & "Maximum predicted energy drift to reduce the timestep",
400 & "/increase split number:",edriftmax
401 write (iout,'(a60,f10.5)')
402 & "Maximum velocity threshold to reduce velocities:",dvmax
403 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
404 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
405 if (rattle) write (iout,'(a60)')
406 & "Rattle algorithm used to constrain the virtual bonds"
410 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
411 call reada(controlcard,"RWAT",rwat,1.4d0)
412 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
413 surfarea=index(controlcard,"SURFAREA").gt.0
414 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
415 if(me.eq.king.or..not.out1file)then
416 write (iout,'(/a,$)') "Langevin dynamics calculation"
419 & " with direct integration of Langevin equations"
420 else if (lang.eq.2) then
421 write (iout,'(a/)') " with TINKER stochasic MD integrator"
422 else if (lang.eq.3) then
423 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
424 else if (lang.eq.4) then
425 write (iout,'(a/)') " in overdamped mode"
427 write (iout,'(//a,i5)')
428 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
431 write (iout,'(a60,f10.5)') "Temperature:",t_bath
432 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
433 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
434 write (iout,'(a60,f10.5)')
435 & "Scaling factor of the friction forces:",scal_fric
436 if (surfarea) write (iout,'(2a,i10,a)')
437 & "Friction coefficients will be scaled by solvent-accessible",
438 & " surface area every",reset_fricmat," steps."
440 c Calculate friction coefficients and bounds of stochastic forces
441 eta=6*pi*cPoise*etawat
442 if(me.eq.king.or..not.out1file)
443 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
445 gamp=scal_fric*(pstok+rwat)*eta
446 stdfp=dsqrt(2*Rb*t_bath/d_time)
448 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
449 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
451 if(me.eq.king.or..not.out1file)then
452 write (iout,'(/2a/)')
453 & "Radii of site types and friction coefficients and std's of",
454 & " stochastic forces of fully exposed sites"
455 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
457 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
458 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
462 if(me.eq.king.or..not.out1file)then
463 write (iout,'(a)') "Berendsen bath calculation"
464 write (iout,'(a60,f10.5)') "Temperature:",t_bath
465 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
467 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
468 & count_reset_moment," steps"
470 & write (iout,'(a,i10,a)')
471 & "Velocities will be reset at random every",count_reset_vel,
475 if(me.eq.king.or..not.out1file)
476 & write (iout,'(a31)') "Microcanonical mode calculation"
478 if(me.eq.king.or..not.out1file)then
479 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
481 write(iout,*) "MD running with constraints."
482 write(iout,*) "Equilibration time ", eq_time, " mtus."
483 write(iout,*) "Constraining ", nfrag," fragments."
484 write(iout,*) "Length of each fragment, weight and q0:"
486 write (iout,*) "Set of restraints #",iset
488 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
489 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
491 write(iout,*) "constraints between ", npair, "fragments."
492 write(iout,*) "constraint pairs, weights and q0:"
494 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
495 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
497 write(iout,*) "angle constraints within ", nfrag_back,
498 & "backbone fragments."
499 write(iout,*) "fragment, weights:"
501 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
502 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
503 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
506 iset=mod(kolor,nset)+1
509 if(me.eq.king.or..not.out1file)
510 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
513 c------------------------------------------------------------------------------
516 C Read molecular data.
518 implicit real*8 (a-h,o-z)
524 include 'COMMON.IOUNITS'
527 include 'COMMON.INTERACT'
528 include 'COMMON.LOCAL'
529 include 'COMMON.NAMES'
530 include 'COMMON.CHAIN'
531 include 'COMMON.FFIELD'
532 include 'COMMON.SBRIDGE'
533 include 'COMMON.HEADER'
534 include 'COMMON.CONTROL'
535 include 'COMMON.DBASE'
536 include 'COMMON.THREAD'
537 include 'COMMON.CONTACTS'
538 include 'COMMON.TORCNSTR'
539 include 'COMMON.TIME1'
540 include 'COMMON.BOUNDS'
542 include 'COMMON.SETUP'
543 character*4 sequence(maxres)
545 double precision x(maxvar)
546 character*256 pdbfile
547 character*320 weightcard
548 character*80 weightcard_t,ucase
549 dimension itype_pdb(maxres)
550 common /pizda/ itype_pdb
551 logical seq_comp,fail
552 double precision energia(0:n_ene)
558 C Read weights of the subsequent energy terms.
559 call card_concat(weightcard)
560 call reada(weightcard,'WLONG',wlong,1.0D0)
561 call reada(weightcard,'WSC',wsc,wlong)
562 call reada(weightcard,'WSCP',wscp,wlong)
563 call reada(weightcard,'WELEC',welec,1.0D0)
564 call reada(weightcard,'WVDWPP',wvdwpp,welec)
565 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
566 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
567 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
568 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
569 call reada(weightcard,'WTURN3',wturn3,1.0D0)
570 call reada(weightcard,'WTURN4',wturn4,1.0D0)
571 call reada(weightcard,'WTURN6',wturn6,1.0D0)
572 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
573 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
574 call reada(weightcard,'WBOND',wbond,1.0D0)
575 call reada(weightcard,'WTOR',wtor,1.0D0)
576 call reada(weightcard,'WTORD',wtor_d,1.0D0)
577 call reada(weightcard,'WANG',wang,1.0D0)
578 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
579 call reada(weightcard,'SCAL14',scal14,0.4D0)
580 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
581 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
582 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
583 call reada(weightcard,'TEMP0',temp0,300.0d0)
584 if (index(weightcard,'SOFT').gt.0) ipot=6
585 C 12/1/95 Added weight for the multi-body term WCORR
586 call reada(weightcard,'WCORRH',wcorr,1.0D0)
587 if (wcorr4.gt.0.0d0) wcorr=wcorr4
607 if(me.eq.king.or..not.out1file)
608 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
609 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
611 10 format (/'Energy-term weights (unscaled):'//
612 & 'WSCC= ',f10.6,' (SC-SC)'/
613 & 'WSCP= ',f10.6,' (SC-p)'/
614 & 'WELEC= ',f10.6,' (p-p electr)'/
615 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
616 & 'WBOND= ',f10.6,' (stretching)'/
617 & 'WANG= ',f10.6,' (bending)'/
618 & 'WSCLOC= ',f10.6,' (SC local)'/
619 & 'WTOR= ',f10.6,' (torsional)'/
620 & 'WTORD= ',f10.6,' (double torsional)'/
621 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
622 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
623 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
624 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
625 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
626 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
627 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
628 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
629 & 'WTURN6= ',f10.6,' (turns, 6th order)')
630 if(me.eq.king.or..not.out1file)then
631 if (wcorr4.gt.0.0d0) then
632 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
633 & 'between contact pairs of peptide groups'
634 write (iout,'(2(a,f5.3/))')
635 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
636 & 'Range of quenching the correlation terms:',2*delt_corr
637 else if (wcorr.gt.0.0d0) then
638 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
639 & 'between contact pairs of peptide groups'
641 write (iout,'(a,f8.3)')
642 & 'Scaling factor of 1,4 SC-p interactions:',scal14
643 write (iout,'(a,f8.3)')
644 & 'General scaling factor of SC-p interactions:',scalscp
646 r0_corr=cutoff_corr-delt_corr
648 aad(i,1)=scalscp*aad(i,1)
649 aad(i,2)=scalscp*aad(i,2)
650 bad(i,1)=scalscp*bad(i,1)
651 bad(i,2)=scalscp*bad(i,2)
653 call rescale_weights(t_bath)
654 if(me.eq.king.or..not.out1file)
655 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
656 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
658 22 format (/'Energy-term weights (scaled):'//
659 & 'WSCC= ',f10.6,' (SC-SC)'/
660 & 'WSCP= ',f10.6,' (SC-p)'/
661 & 'WELEC= ',f10.6,' (p-p electr)'/
662 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
663 & 'WBOND= ',f10.6,' (stretching)'/
664 & 'WANG= ',f10.6,' (bending)'/
665 & 'WSCLOC= ',f10.6,' (SC local)'/
666 & 'WTOR= ',f10.6,' (torsional)'/
667 & 'WTORD= ',f10.6,' (double torsional)'/
668 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
669 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
670 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
671 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
672 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
673 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
674 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
675 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
676 & 'WTURN6= ',f10.6,' (turns, 6th order)')
677 if(me.eq.king.or..not.out1file)
678 & write (iout,*) "Reference temperature for weights calculation:",
680 call reada(weightcard,"D0CM",d0cm,3.78d0)
681 call reada(weightcard,"AKCM",akcm,15.1d0)
682 call reada(weightcard,"AKTH",akth,11.0d0)
683 call reada(weightcard,"AKCT",akct,12.0d0)
684 call reada(weightcard,"V1SS",v1ss,-1.08d0)
685 call reada(weightcard,"V2SS",v2ss,7.61d0)
686 call reada(weightcard,"V3SS",v3ss,13.7d0)
687 call reada(weightcard,"EBR",ebr,-5.50D0)
688 if(me.eq.king.or..not.out1file) then
689 write (iout,*) "Parameters of the SS-bond potential:"
690 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
692 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
693 write (iout,*) "EBR",ebr
694 c print *,'indpdb=',indpdb,' pdbref=',pdbref
696 if (indpdb.gt.0 .or. pdbref) then
697 read(inp,'(a)') pdbfile
698 if(me.eq.king.or..not.out1file)
699 & write (iout,'(2a)') 'PDB data will be read from file ',
700 & pdbfile(:ilen(pdbfile))
701 open(ipdbin,file=pdbfile,status='old',err=33)
703 33 write (iout,'(a)') 'Error opening PDB file.'
706 c print *,'Begin reading pdb data'
708 c print *,'Finished reading pdb data'
709 if(me.eq.king.or..not.out1file)
710 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
711 & ' nstart_sup=',nstart_sup
713 itype_pdb(i)=itype(i)
717 nct=nstart_sup+nsup-1
718 call contact(.false.,ncont_ref,icont_ref,co)
721 if(me.eq.king.or..not.out1file)
722 & write(iout,*)'Adding sidechains'
726 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
729 do while (fail.and.nsi.le.maxsi)
730 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
733 if(fail) write(iout,*)'Adding sidechain failed for res ',
734 & i,' after ',nsi,' trials'
739 if (indpdb.eq.0) then
740 C Read sequence if not taken from the pdb file.
742 c print *,'nres=',nres
743 if (iscode.gt.0) then
744 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
746 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
748 C Convert sequence to numeric code
750 itype(i)=rescode(i,sequence(i),iscode)
752 C Assign initial virtual bond lengths
758 vbld(i+nres)=dsc(iabs(itype(i)))
759 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
760 c write (iout,*) "i",i," itype",itype(i),
761 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
765 c print '(20i4)',(itype(i),i=1,nres)
768 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
770 if (itype(i).eq.ntyp1) then
774 else if (iabs(itype(i+1)).ne.20) then
776 else if (iabs(itype(i)).ne.20) then
783 if(me.eq.king.or..not.out1file)then
784 write (iout,*) "ITEL"
786 write (iout,*) i,itype(i),itel(i)
788 print *,'Call Read_Bridge.'
791 C 8/13/98 Set limits to generating the dihedral angles
796 read (inp,*) ndih_constr
797 if (ndih_constr.gt.0) then
799 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
800 if(me.eq.king.or..not.out1file)then
802 & 'There are',ndih_constr,' constraints on phi angles.'
804 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
808 phi0(i)=deg2rad*phi0(i)
809 drange(i)=deg2rad*drange(i)
811 if(me.eq.king.or..not.out1file)
812 & write (iout,*) 'FTORS',ftors
815 phibound(1,ii) = phi0(i)-drange(i)
816 phibound(2,ii) = phi0(i)+drange(i)
823 write (iout,'(a)') 'Boundaries in phi angle sampling:'
825 write (iout,'(a3,i5,2f10.1)')
826 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
832 cd print *,'NNT=',NNT,' NCT=',NCT
833 if (itype(1).eq.ntyp1) nnt=2
834 if (itype(nres).eq.ntyp1) nct=nct-1
836 if(me.eq.king.or..not.out1file)
837 & write (iout,'(a,i3)') 'nsup=',nsup
839 if (nsup.le.(nct-nnt+1)) then
840 do i=0,nct-nnt+1-nsup
841 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
847 & 'Error - sequences to be superposed do not match.'
850 do i=0,nsup-(nct-nnt+1)
851 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
853 nstart_sup=nstart_sup+i
859 & 'Error - sequences to be superposed do not match.'
862 if (nsup.eq.0) nsup=nct-nnt
863 if (nstart_sup.eq.0) nstart_sup=nnt
864 if (nstart_seq.eq.0) nstart_seq=nnt
865 if(me.eq.king.or..not.out1file)
866 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
867 & ' nstart_seq=',nstart_seq
869 c--- Zscore rms -------
870 if (nz_start.eq.0) nz_start=nnt
871 if (nz_end.eq.0 .and. nsup.gt.0) then
873 else if (nz_end.eq.0) then
876 if(me.eq.king.or..not.out1file)then
877 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
878 write (iout,*) 'IZ_SC=',iz_sc
880 c----------------------
883 if (.not.pdbref) then
884 call read_angles(inp,*38)
886 38 write (iout,'(a)') 'Error reading reference structure.'
888 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
889 stop 'Error reading reference structure'
893 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
903 call contact(.true.,ncont_ref,icont_ref,co)
905 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
907 if (constr_dist.gt.0) call read_dist_constr
908 write (iout,*) "After read_dist_constr nhpb",nhpb
910 if(me.eq.king.or..not.out1file)
911 & write (iout,*) 'Contact order:',co
913 if(me.eq.king.or..not.out1file)
914 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
917 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
919 if(me.eq.king.or..not.out1file)
920 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
921 & icont_ref(1,i),' ',
922 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
926 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
927 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
928 & modecalc.ne.10) then
929 C If input structure hasn't been supplied from the PDB file read or generate
931 if (iranconf.eq.0 .and. .not. extconf) then
932 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
933 & write (iout,'(a)') 'Initial geometry will be read in.'
935 read(inp,'(8f10.5)',end=36,err=36)
936 & ((c(l,k),l=1,3),k=1,nres),
937 & ((c(l,k+nres),l=1,3),k=nnt,nct)
938 write (iout,*) "Exit READ_CART"
939 write (iout,'(8f10.5)')
940 & ((c(l,k),l=1,3),k=1,nres),
941 & ((c(l,k+nres),l=1,3),k=nnt,nct)
942 call int_from_cart1(.true.)
943 write (iout,*) "Finish INT_TO_CART"
946 dc(j,i)=c(j,i+1)-c(j,i)
947 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
951 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
953 dc(j,i+nres)=c(j,i+nres)-c(j,i)
954 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
960 call read_angles(inp,*36)
963 36 write (iout,'(a)') 'Error reading angle file.'
965 call mpi_finalize( MPI_COMM_WORLD,IERR )
967 stop 'Error reading angle file.'
969 else if (extconf) then
970 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
971 & write (iout,'(a)') 'Extended chain initial geometry.'
973 theta(i)=90d0*deg2rad
979 alph(i)=110d0*deg2rad
982 omeg(i)=-120d0*deg2rad
983 if (itype(i).le.0) omeg(i)=-omeg(i)
986 if(me.eq.king.or..not.out1file)
987 & write (iout,'(a)') 'Random-generated initial geometry.'
991 if (me.eq.king .or. fg_rank.eq.0 .and. (
992 & modecalc.eq.12 .or. modecalc.eq.14) ) then
996 call gen_rand_conf(itmp,*30)
998 30 write (iout,*) 'Failed to generate random conformation',
1000 write (*,*) 'Processor:',me,
1001 & ' Failed to generate random conformation',
1011 write (iout,'(a,i3,a)') 'Processor:',me,
1012 & ' error in generating random conformation.'
1013 write (*,'(a,i3,a)') 'Processor:',me,
1014 & ' error in generating random conformation.'
1017 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1023 call gen_rand_conf(itmp,*30)
1025 30 write (iout,*) 'Failed to generate random conformation',
1026 & ', itrial=',itrial
1027 write (*,*) 'Failed to generate random conformation',
1028 & ', itrial=',itrial
1030 write (iout,'(a,i3,a)') 'Processor:',me,
1031 & ' error in generating random conformation.'
1032 write (*,'(a,i3,a)') 'Processor:',me,
1033 & ' error in generating random conformation.'
1038 elseif (modecalc.eq.4) then
1039 read (inp,'(a)') intinname
1040 open (intin,file=intinname,status='old',err=333)
1041 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1042 & write (iout,'(a)') 'intinname',intinname
1043 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1045 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1047 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1049 stop 'Error opening angle file.'
1053 C Generate distance constraints, if the PDB structure is to be regularized.
1054 if (nthread.gt.0) then
1055 call read_threadbase
1058 if (me.eq.king .or. .not. out1file)
1060 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1061 write (iout,'(/a,i3,a)')
1062 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1063 write (iout,'(20i4)') (iss(i),i=1,ns)
1064 write (iout,'(/a/)') 'Pre-formed links are:'
1070 if (me.eq.king.or..not.out1file)
1071 & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1072 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1077 if (i2ndstr.gt.0) call secstrp2dihc
1078 c call geom_to_var(nvar,x)
1079 c call etotal(energia(0))
1080 c call enerprint(energia(0))
1081 c call briefout(0,etot)
1083 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1084 cd write (iout,'(a)') 'Variable list:'
1085 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1087 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1088 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1089 & 'Processor',myrank,': end reading molecular data.'
1093 c--------------------------------------------------------------------------
1094 logical function seq_comp(itypea,itypeb,length)
1096 integer length,itypea(length),itypeb(length)
1099 if (itypea(i).ne.itypeb(i)) then
1107 c-----------------------------------------------------------------------------
1108 subroutine read_bridge
1109 C Read information about disulfide bridges.
1110 implicit real*8 (a-h,o-z)
1111 include 'DIMENSIONS'
1115 include 'COMMON.IOUNITS'
1116 include 'COMMON.GEO'
1117 include 'COMMON.VAR'
1118 include 'COMMON.INTERACT'
1119 include 'COMMON.LOCAL'
1120 include 'COMMON.NAMES'
1121 include 'COMMON.CHAIN'
1122 include 'COMMON.FFIELD'
1123 include 'COMMON.SBRIDGE'
1124 include 'COMMON.HEADER'
1125 include 'COMMON.CONTROL'
1126 include 'COMMON.DBASE'
1127 include 'COMMON.THREAD'
1128 include 'COMMON.TIME1'
1129 include 'COMMON.SETUP'
1130 C Read bridging residues.
1131 read (inp,*) ns,(iss(i),i=1,ns)
1133 if(me.eq.king.or..not.out1file)
1134 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1135 C Check whether the specified bridging residues are cystines.
1137 if (itype(iss(i)).ne.1) then
1138 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1139 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
1140 & ' can form a disulfide bridge?!!!'
1141 write (*,'(2a,i3,a)')
1142 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
1143 & ' can form a disulfide bridge?!!!'
1145 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1150 C Read preformed bridges.
1152 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1153 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1156 C Check if the residues involved in bridges are in the specified list of
1157 C bridging residues.
1160 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1161 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1162 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1163 & ' contains residues present in other pairs.'
1164 write (*,'(a,i3,a)') 'Disulfide pair',i,
1165 & ' contains residues present in other pairs.'
1167 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1173 if (ihpb(i).eq.iss(j)) goto 10
1175 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1178 if (jhpb(i).eq.iss(j)) goto 20
1180 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1186 ihpb(i)=ihpb(i)+nres
1187 jhpb(i)=jhpb(i)+nres
1193 c----------------------------------------------------------------------------
1194 subroutine read_x(kanal,*)
1195 implicit real*8 (a-h,o-z)
1196 include 'DIMENSIONS'
1197 include 'COMMON.GEO'
1198 include 'COMMON.VAR'
1199 include 'COMMON.CHAIN'
1200 include 'COMMON.IOUNITS'
1201 include 'COMMON.CONTROL'
1202 include 'COMMON.LOCAL'
1203 include 'COMMON.INTERACT'
1204 c Read coordinates from input
1206 read(kanal,'(8f10.5)',end=10,err=10)
1207 & ((c(l,k),l=1,3),k=1,nres),
1208 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1211 c(j,2*nres)=c(j,nres)
1213 call int_from_cart1(.false.)
1216 dc(j,i)=c(j,i+1)-c(j,i)
1217 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1221 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1223 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1224 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1232 c----------------------------------------------------------------------------
1233 subroutine read_threadbase
1234 implicit real*8 (a-h,o-z)
1235 include 'DIMENSIONS'
1236 include 'COMMON.IOUNITS'
1237 include 'COMMON.GEO'
1238 include 'COMMON.VAR'
1239 include 'COMMON.INTERACT'
1240 include 'COMMON.LOCAL'
1241 include 'COMMON.NAMES'
1242 include 'COMMON.CHAIN'
1243 include 'COMMON.FFIELD'
1244 include 'COMMON.SBRIDGE'
1245 include 'COMMON.HEADER'
1246 include 'COMMON.CONTROL'
1247 include 'COMMON.DBASE'
1248 include 'COMMON.THREAD'
1249 include 'COMMON.TIME1'
1250 C Read pattern database for threading.
1251 read (icbase,*) nseq
1253 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1254 & nres_base(2,i),nres_base(3,i)
1255 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1257 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1258 c & nres_base(2,i),nres_base(3,i)
1259 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1263 if (weidis.eq.0.0D0) weidis=0.1D0
1272 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1273 write (iout,'(a,i5)') 'nexcl: ',nexcl
1274 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1277 c------------------------------------------------------------------------------
1278 subroutine setup_var
1279 implicit real*8 (a-h,o-z)
1280 include 'DIMENSIONS'
1281 include 'COMMON.IOUNITS'
1282 include 'COMMON.GEO'
1283 include 'COMMON.VAR'
1284 include 'COMMON.INTERACT'
1285 include 'COMMON.LOCAL'
1286 include 'COMMON.NAMES'
1287 include 'COMMON.CHAIN'
1288 include 'COMMON.FFIELD'
1289 include 'COMMON.SBRIDGE'
1290 include 'COMMON.HEADER'
1291 include 'COMMON.CONTROL'
1292 include 'COMMON.DBASE'
1293 include 'COMMON.THREAD'
1294 include 'COMMON.TIME1'
1295 C Set up variable list.
1301 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1303 ialph(i,1)=nvar+nside
1307 if (indphi.gt.0) then
1309 else if (indback.gt.0) then
1314 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1317 c----------------------------------------------------------------------------
1318 subroutine gen_dist_constr
1319 C Generate CA distance constraints.
1320 implicit real*8 (a-h,o-z)
1321 include 'DIMENSIONS'
1322 include 'COMMON.IOUNITS'
1323 include 'COMMON.GEO'
1324 include 'COMMON.VAR'
1325 include 'COMMON.INTERACT'
1326 include 'COMMON.LOCAL'
1327 include 'COMMON.NAMES'
1328 include 'COMMON.CHAIN'
1329 include 'COMMON.FFIELD'
1330 include 'COMMON.SBRIDGE'
1331 include 'COMMON.HEADER'
1332 include 'COMMON.CONTROL'
1333 include 'COMMON.DBASE'
1334 include 'COMMON.THREAD'
1335 include 'COMMON.TIME1'
1336 dimension itype_pdb(maxres)
1337 common /pizda/ itype_pdb
1339 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1340 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1341 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1343 do i=nstart_sup,nstart_sup+nsup-1
1344 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1345 cd & ' seq_pdb', restyp(itype_pdb(i))
1346 do j=i+2,nstart_sup+nsup-1
1348 ihpb(nhpb)=i+nstart_seq-nstart_sup
1349 jhpb(nhpb)=j+nstart_seq-nstart_sup
1351 dhpb(nhpb)=dist(i,j)
1354 cd write (iout,'(a)') 'Distance constraints:'
1359 cd if (ii.gt.nres) then
1364 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1365 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1366 cd & dhpb(i),forcon(i)
1370 c----------------------------------------------------------------------------
1372 implicit real*8 (a-h,o-z)
1373 include 'DIMENSIONS'
1374 include 'COMMON.MAP'
1375 include 'COMMON.IOUNITS'
1376 character*3 angid(4) /'THE','PHI','ALP','OME'/
1377 character*80 mapcard,ucase
1379 read (inp,'(a)') mapcard
1380 mapcard=ucase(mapcard)
1381 if (index(mapcard,'PHI').gt.0) then
1383 else if (index(mapcard,'THE').gt.0) then
1385 else if (index(mapcard,'ALP').gt.0) then
1387 else if (index(mapcard,'OME').gt.0) then
1390 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1391 stop 'Error - illegal variable spec in MAP card.'
1393 call readi (mapcard,'RES1',res1(imap),0)
1394 call readi (mapcard,'RES2',res2(imap),0)
1395 if (res1(imap).eq.0) then
1396 res1(imap)=res2(imap)
1397 else if (res2(imap).eq.0) then
1398 res2(imap)=res1(imap)
1400 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1402 & 'Error - illegal definition of variable group in MAP.'
1403 stop 'Error - illegal definition of variable group in MAP.'
1405 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1406 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1407 call readi(mapcard,'NSTEP',nstep(imap),0)
1408 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1410 & 'Illegal boundary and/or step size specification in MAP.'
1411 stop 'Illegal boundary and/or step size specification in MAP.'
1416 c----------------------------------------------------------------------------
1418 implicit real*8 (a-h,o-z)
1419 include 'DIMENSIONS'
1420 include 'COMMON.IOUNITS'
1421 include 'COMMON.GEO'
1422 include 'COMMON.CSA'
1423 include 'COMMON.BANK'
1424 include 'COMMON.CONTROL'
1426 character*620 mcmcard
1427 call card_concat(mcmcard)
1429 call readi(mcmcard,'NCONF',nconf,50)
1430 call readi(mcmcard,'NADD',nadd,0)
1431 call readi(mcmcard,'JSTART',jstart,1)
1432 call readi(mcmcard,'JEND',jend,1)
1433 call readi(mcmcard,'NSTMAX',nstmax,500000)
1434 call readi(mcmcard,'N0',n0,1)
1435 call readi(mcmcard,'N1',n1,6)
1436 call readi(mcmcard,'N2',n2,4)
1437 call readi(mcmcard,'N3',n3,0)
1438 call readi(mcmcard,'N4',n4,0)
1439 call readi(mcmcard,'N5',n5,0)
1440 call readi(mcmcard,'N6',n6,10)
1441 call readi(mcmcard,'N7',n7,0)
1442 call readi(mcmcard,'N8',n8,0)
1443 call readi(mcmcard,'N9',n9,0)
1444 call readi(mcmcard,'N14',n14,0)
1445 call readi(mcmcard,'N15',n15,0)
1446 call readi(mcmcard,'N16',n16,0)
1447 call readi(mcmcard,'N17',n17,0)
1448 call readi(mcmcard,'N18',n18,0)
1450 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1452 call readi(mcmcard,'NDIFF',ndiff,2)
1453 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1454 call readi(mcmcard,'IS1',is1,1)
1455 call readi(mcmcard,'IS2',is2,8)
1456 call readi(mcmcard,'NRAN0',nran0,4)
1457 call readi(mcmcard,'NRAN1',nran1,2)
1458 call readi(mcmcard,'IRR',irr,1)
1459 call readi(mcmcard,'NSEED',nseed,20)
1460 call readi(mcmcard,'NTOTAL',ntotal,10000)
1461 call reada(mcmcard,'CUT1',cut1,2.0d0)
1462 call reada(mcmcard,'CUT2',cut2,5.0d0)
1463 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1464 call readi(mcmcard,'ICMAX',icmax,3)
1465 call readi(mcmcard,'IRESTART',irestart,0)
1466 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1469 call reada(mcmcard,'DELE',dele,20.0d0)
1470 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1471 call readi(mcmcard,'IREF',iref,0)
1472 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1473 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1474 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1475 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1476 write (iout,*) "NCONF_IN",nconf_in
1479 c----------------------------------------------------------------------------
1480 cfmc subroutine mcmfread
1481 cfmc implicit real*8 (a-h,o-z)
1482 cfmc include 'DIMENSIONS'
1483 cfmc include 'COMMON.MCMF'
1484 cfmc include 'COMMON.IOUNITS'
1485 cfmc include 'COMMON.GEO'
1486 cfmc character*80 ucase
1487 cfmc character*620 mcmcard
1488 cfmc call card_concat(mcmcard)
1490 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1491 cfmc write(iout,*)'MAXRANT=',maxrant
1492 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1493 cfmc write(iout,*)'MAXFAM=',maxfam
1494 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1495 cfmc write(iout,*)'NNET1=',nnet1
1496 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1497 cfmc write(iout,*)'NNET2=',nnet2
1498 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1499 cfmc write(iout,*)'NNET3=',nnet3
1500 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1501 cfmc write(iout,*)'ILASTT=',ilastt
1502 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1503 cfmc write(iout,*)'MAXSTR=',maxstr
1504 cfmc maxstr_f=maxstr/maxfam
1505 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1506 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1507 cfmc write(iout,*)'NMCMF=',nmcmf
1508 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1509 cfmc write(iout,*)'IFOCUS=',ifocus
1510 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1511 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1512 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1513 cfmc write(iout,*)'INTPRT=',intprt
1514 cfmc call readi(mcmcard,'IPRT',iprt,100)
1515 cfmc write(iout,*)'IPRT=',iprt
1516 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1517 cfmc write(iout,*)'IMAXTR=',imaxtr
1518 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1519 cfmc write(iout,*)'MAXEVEN=',maxeven
1520 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1521 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1522 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1523 cfmc write(iout,*)'INIMIN=',inimin
1524 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1525 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1526 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1527 cfmc write(iout,*)'NTHREAD=',nthread
1528 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1529 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1530 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1531 cfmc write(iout,*)'MAXPERT=',maxpert
1532 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1533 cfmc write(iout,*)'IRMSD=',irmsd
1534 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1535 cfmc write(iout,*)'DENEMIN=',denemin
1536 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1537 cfmc write(iout,*)'RCUT1S=',rcut1s
1538 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1539 cfmc write(iout,*)'RCUT1E=',rcut1e
1540 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1541 cfmc write(iout,*)'RCUT2S=',rcut2s
1542 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1543 cfmc write(iout,*)'RCUT2E=',rcut2e
1544 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1545 cfmc write(iout,*)'DPERT1=',d_pert1
1546 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1547 cfmc write(iout,*)'DPERT1A=',d_pert1a
1548 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1549 cfmc write(iout,*)'DPERT2=',d_pert2
1550 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1551 cfmc write(iout,*)'DPERT2A=',d_pert2a
1552 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1553 cfmc write(iout,*)'DPERT2B=',d_pert2b
1554 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1555 cfmc write(iout,*)'DPERT2C=',d_pert2c
1556 cfmc d_pert1=deg2rad*d_pert1
1557 cfmc d_pert1a=deg2rad*d_pert1a
1558 cfmc d_pert2=deg2rad*d_pert2
1559 cfmc d_pert2a=deg2rad*d_pert2a
1560 cfmc d_pert2b=deg2rad*d_pert2b
1561 cfmc d_pert2c=deg2rad*d_pert2c
1562 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1563 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1564 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1565 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1566 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1567 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1568 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1569 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1570 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1571 cfmc write(iout,*)'RCUTINI=',rcutini
1572 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1573 cfmc write(iout,*)'GRAT=',grat
1574 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1575 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1579 c----------------------------------------------------------------------------
1581 implicit real*8 (a-h,o-z)
1582 include 'DIMENSIONS'
1583 include 'COMMON.MCM'
1584 include 'COMMON.MCE'
1585 include 'COMMON.IOUNITS'
1587 character*320 mcmcard
1588 call card_concat(mcmcard)
1589 call readi(mcmcard,'MAXACC',maxacc,100)
1590 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1591 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1592 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1593 call readi(mcmcard,'MAXREPM',maxrepm,200)
1594 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1595 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1596 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1597 call reada(mcmcard,'E_UP',e_up,5.0D0)
1598 call reada(mcmcard,'DELTE',delte,0.1D0)
1599 call readi(mcmcard,'NSWEEP',nsweep,5)
1600 call readi(mcmcard,'NSTEPH',nsteph,0)
1601 call readi(mcmcard,'NSTEPC',nstepc,0)
1602 call reada(mcmcard,'TMIN',tmin,298.0D0)
1603 call reada(mcmcard,'TMAX',tmax,298.0D0)
1604 call readi(mcmcard,'NWINDOW',nwindow,0)
1605 call readi(mcmcard,'PRINT_MC',print_mc,0)
1606 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1607 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1608 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1609 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1610 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1611 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1612 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1613 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1614 if (nwindow.gt.0) then
1615 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1617 winlen(i)=winend(i)-winstart(i)+1
1620 if (tmax.lt.tmin) tmax=tmin
1621 if (tmax.eq.tmin) then
1625 if (nstepc.gt.0 .and. nsteph.gt.0) then
1626 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1627 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1629 C Probabilities of different move types
1630 sumpro_type(0)=0.0D0
1631 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1632 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1633 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1634 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1635 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1636 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1637 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1639 print *,'i',i,' sumprotype',sumpro_type(i)
1640 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1641 print *,'i',i,' sumprotype',sumpro_type(i)
1645 c----------------------------------------------------------------------------
1646 subroutine read_minim
1647 implicit real*8 (a-h,o-z)
1648 include 'DIMENSIONS'
1649 include 'COMMON.MINIM'
1650 include 'COMMON.IOUNITS'
1652 character*320 minimcard
1653 call card_concat(minimcard)
1654 call readi(minimcard,'MAXMIN',maxmin,2000)
1655 call readi(minimcard,'MAXFUN',maxfun,5000)
1656 call readi(minimcard,'MINMIN',minmin,maxmin)
1657 call readi(minimcard,'MINFUN',minfun,maxmin)
1658 call reada(minimcard,'TOLF',tolf,1.0D-2)
1659 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1660 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1661 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1662 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1663 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1664 & 'Options in energy minimization:'
1665 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1666 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1667 & 'MinMin:',MinMin,' MinFun:',MinFun,
1668 & ' TolF:',TolF,' RTolF:',RTolF
1671 c----------------------------------------------------------------------------
1672 subroutine read_angles(kanal,*)
1673 implicit real*8 (a-h,o-z)
1674 include 'DIMENSIONS'
1675 include 'COMMON.GEO'
1676 include 'COMMON.VAR'
1677 include 'COMMON.CHAIN'
1678 include 'COMMON.IOUNITS'
1679 include 'COMMON.CONTROL'
1680 c Read angles from input
1682 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1683 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1684 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1685 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1688 c 9/7/01 avoid 180 deg valence angle
1689 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1691 theta(i)=deg2rad*theta(i)
1692 phi(i)=deg2rad*phi(i)
1693 alph(i)=deg2rad*alph(i)
1694 omeg(i)=deg2rad*omeg(i)
1699 c----------------------------------------------------------------------------
1700 subroutine reada(rekord,lancuch,wartosc,default)
1702 character*(*) rekord,lancuch
1703 double precision wartosc,default
1706 iread=index(rekord,lancuch)
1707 if (iread.eq.0) then
1711 iread=iread+ilen(lancuch)+1
1712 read (rekord(iread:),*,err=10,end=10) wartosc
1717 c----------------------------------------------------------------------------
1718 subroutine readi(rekord,lancuch,wartosc,default)
1720 character*(*) rekord,lancuch
1721 integer wartosc,default
1724 iread=index(rekord,lancuch)
1725 if (iread.eq.0) then
1729 iread=iread+ilen(lancuch)+1
1730 read (rekord(iread:),*,err=10,end=10) wartosc
1735 c----------------------------------------------------------------------------
1736 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1739 integer tablica(dim),default
1740 character*(*) rekord,lancuch
1747 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1748 if (iread.eq.0) return
1749 iread=iread+ilen(lancuch)+1
1750 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1753 c----------------------------------------------------------------------------
1754 subroutine multreada(rekord,lancuch,tablica,dim,default)
1757 double precision tablica(dim),default
1758 character*(*) rekord,lancuch
1765 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1766 if (iread.eq.0) return
1767 iread=iread+ilen(lancuch)+1
1768 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1771 c----------------------------------------------------------------------------
1772 subroutine openunits
1773 implicit real*8 (a-h,o-z)
1774 include 'DIMENSIONS'
1777 character*16 form,nodename
1780 include 'COMMON.SETUP'
1781 include 'COMMON.IOUNITS'
1783 include 'COMMON.CONTROL'
1784 integer lenpre,lenpot,ilen,lentmp
1786 character*3 out1file_text,ucase
1789 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1790 call getenv_loc("PREFIX",prefix)
1792 call getenv_loc("POT",pot)
1793 call getenv_loc("DIRTMP",tmpdir)
1794 call getenv_loc("CURDIR",curdir)
1795 call getenv_loc("OUT1FILE",out1file_text)
1796 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1797 out1file_text=ucase(out1file_text)
1798 if (out1file_text(1:1).eq."Y") then
1801 out1file=fg_rank.gt.0
1806 if (lentmp.gt.0) then
1807 write (*,'(80(1h!))')
1808 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1809 write (*,'(80(1h!))')
1810 write (*,*)"All output files will be on node /tmp directory."
1812 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1813 if (me.eq.king) then
1814 write (*,*) "The master node is ",nodename
1815 else if (fg_rank.eq.0) then
1816 write (*,*) "I am the CG slave node ",nodename
1818 write (*,*) "I am the FG slave node ",nodename
1821 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1822 lenpre = lentmp+lenpre+1
1824 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1825 C Get the names and open the input files
1826 #if defined(WINIFL) || defined(WINPGI)
1827 open(1,file=pref_orig(:ilen(pref_orig))//
1828 & '.inp',status='old',readonly,shared)
1829 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1830 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1831 C Get parameter filenames and open the parameter files.
1832 call getenv_loc('BONDPAR',bondname)
1833 open (ibond,file=bondname,status='old',readonly,shared)
1834 call getenv_loc('THETPAR',thetname)
1835 open (ithep,file=thetname,status='old',readonly,shared)
1836 call getenv_loc('ROTPAR',rotname)
1837 open (irotam,file=rotname,status='old',readonly,shared)
1838 call getenv_loc('TORPAR',torname)
1839 open (itorp,file=torname,status='old',readonly,shared)
1840 call getenv_loc('TORDPAR',tordname)
1841 open (itordp,file=tordname,status='old',readonly,shared)
1842 call getenv_loc('FOURIER',fouriername)
1843 open (ifourier,file=fouriername,status='old',readonly,shared)
1844 call getenv_loc('ELEPAR',elename)
1845 open (ielep,file=elename,status='old',readonly,shared)
1846 call getenv_loc('SIDEPAR',sidename)
1847 open (isidep,file=sidename,status='old',readonly,shared)
1848 #elif (defined CRAY) || (defined AIX)
1849 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1851 c print *,"Processor",myrank," opened file 1"
1852 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1853 c print *,"Processor",myrank," opened file 9"
1854 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1855 C Get parameter filenames and open the parameter files.
1856 call getenv_loc('BONDPAR',bondname)
1857 open (ibond,file=bondname,status='old',action='read')
1858 c print *,"Processor",myrank," opened file IBOND"
1859 call getenv_loc('THETPAR',thetname)
1860 open (ithep,file=thetname,status='old',action='read')
1861 c print *,"Processor",myrank," opened file ITHEP"
1862 call getenv_loc('ROTPAR',rotname)
1863 open (irotam,file=rotname,status='old',action='read')
1864 c print *,"Processor",myrank," opened file IROTAM"
1865 call getenv_loc('TORPAR',torname)
1866 open (itorp,file=torname,status='old',action='read')
1867 c print *,"Processor",myrank," opened file ITORP"
1868 call getenv_loc('TORDPAR',tordname)
1869 open (itordp,file=tordname,status='old',action='read')
1870 c print *,"Processor",myrank," opened file ITORDP"
1871 call getenv_loc('SCCORPAR',sccorname)
1872 open (isccor,file=sccorname,status='old',action='read')
1873 c print *,"Processor",myrank," opened file ISCCOR"
1874 call getenv_loc('FOURIER',fouriername)
1875 open (ifourier,file=fouriername,status='old',action='read')
1876 c print *,"Processor",myrank," opened file IFOURIER"
1877 call getenv_loc('ELEPAR',elename)
1878 open (ielep,file=elename,status='old',action='read')
1879 c print *,"Processor",myrank," opened file IELEP"
1880 call getenv_loc('SIDEPAR',sidename)
1881 open (isidep,file=sidename,status='old',action='read')
1882 c print *,"Processor",myrank," opened file ISIDEP"
1883 c print *,"Processor",myrank," opened parameter files"
1885 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1886 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1887 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1888 C Get parameter filenames and open the parameter files.
1889 call getenv_loc('BONDPAR',bondname)
1890 open (ibond,file=bondname,status='old')
1891 call getenv_loc('THETPAR',thetname)
1892 open (ithep,file=thetname,status='old')
1893 call getenv_loc('ROTPAR',rotname)
1894 open (irotam,file=rotname,status='old')
1895 call getenv_loc('TORPAR',torname)
1896 open (itorp,file=torname,status='old')
1897 call getenv_loc('TORDPAR',tordname)
1898 open (itordp,file=tordname,status='old')
1899 call getenv_loc('SCCORPAR',sccorname)
1900 open (isccor,file=sccorname,status='old')
1901 call getenv_loc('FOURIER',fouriername)
1902 open (ifourier,file=fouriername,status='old')
1903 call getenv_loc('ELEPAR',elename)
1904 open (ielep,file=elename,status='old')
1905 call getenv_loc('SIDEPAR',sidename)
1906 open (isidep,file=sidename,status='old')
1908 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1910 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1911 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1912 C Get parameter filenames and open the parameter files.
1913 call getenv_loc('BONDPAR',bondname)
1914 open (ibond,file=bondname,status='old',readonly)
1915 call getenv_loc('THETPAR',thetname)
1916 open (ithep,file=thetname,status='old',readonly)
1917 call getenv_loc('ROTPAR',rotname)
1918 open (irotam,file=rotname,status='old',readonly)
1919 call getenv_loc('TORPAR',torname)
1920 open (itorp,file=torname,status='old',readonly)
1921 call getenv_loc('TORDPAR',tordname)
1922 open (itordp,file=tordname,status='old',readonly)
1923 call getenv_loc('SCCORPAR',sccorname)
1924 open (isccor,file=sccorname,status='old',readonly)
1926 call getenv_loc('THETPARPDB',thetname_pdb)
1927 print *,"thetname_pdb ",thetname_pdb
1928 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1929 print *,ithep_pdb," opened"
1931 call getenv_loc('FOURIER',fouriername)
1932 open (ifourier,file=fouriername,status='old',readonly)
1933 call getenv_loc('ELEPAR',elename)
1934 open (ielep,file=elename,status='old',readonly)
1935 call getenv_loc('SIDEPAR',sidename)
1936 open (isidep,file=sidename,status='old',readonly)
1938 call getenv_loc('ROTPARPDB',rotname_pdb)
1939 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1944 C 8/9/01 In the newest version SCp interaction constants are read from a file
1945 C Use -DOLDSCP to use hard-coded constants instead.
1947 call getenv_loc('SCPPAR',scpname)
1948 #if defined(WINIFL) || defined(WINPGI)
1949 open (iscpp,file=scpname,status='old',readonly,shared)
1950 #elif (defined CRAY) || (defined AIX)
1951 open (iscpp,file=scpname,status='old',action='read')
1953 open (iscpp,file=scpname,status='old')
1955 open (iscpp,file=scpname,status='old',readonly)
1958 call getenv_loc('PATTERN',patname)
1959 #if defined(WINIFL) || defined(WINPGI)
1960 open (icbase,file=patname,status='old',readonly,shared)
1961 #elif (defined CRAY) || (defined AIX)
1962 open (icbase,file=patname,status='old',action='read')
1964 open (icbase,file=patname,status='old')
1966 open (icbase,file=patname,status='old',readonly)
1969 C Open output file only for CG processes
1970 c print *,"Processor",myrank," fg_rank",fg_rank
1971 if (fg_rank.eq.0) then
1973 if (nodes.eq.1) then
1976 npos = dlog10(dfloat(nodes-1))+1
1978 if (npos.lt.3) npos=3
1979 write (liczba,'(i1)') npos
1980 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1982 write (liczba,form) me
1983 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1984 & liczba(:ilen(liczba))
1985 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1987 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1989 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1990 & liczba(:ilen(liczba))//'.mol2'
1991 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1992 & liczba(:ilen(liczba))//'.stat'
1994 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1995 & //liczba(:ilen(liczba))//'.stat')
1996 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1999 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2000 & liczba(:ilen(liczba))//'.const'
2005 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2006 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2007 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2008 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2009 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2011 & call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
2013 rest2name=prefix(:ilen(prefix))//'.rst'
2015 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2018 #if defined(AIX) || defined(PGI)
2019 if (me.eq.king .or. .not. out1file)
2020 & open(iout,file=outname,status='unknown')
2022 if (fg_rank.gt.0) then
2023 write (liczba,'(i3.3)') myrank/nfgtasks
2024 write (ll,'(bz,i3.3)') fg_rank
2025 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2030 open(igeom,file=intname,status='unknown',position='append')
2031 open(ipdb,file=pdbname,status='unknown')
2032 open(imol2,file=mol2name,status='unknown')
2033 open(istat,file=statname,status='unknown',position='append')
2035 c1out open(iout,file=outname,status='unknown')
2038 if (me.eq.king .or. .not.out1file)
2039 & open(iout,file=outname,status='unknown')
2041 if (fg_rank.gt.0) then
2042 write (liczba,'(i3.3)') myrank/nfgtasks
2043 write (ll,'(bz,i3.3)') fg_rank
2044 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2049 open(igeom,file=intname,status='unknown',access='append')
2050 open(ipdb,file=pdbname,status='unknown')
2051 open(imol2,file=mol2name,status='unknown')
2052 open(istat,file=statname,status='unknown',access='append')
2054 c1out open(iout,file=outname,status='unknown')
2057 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2058 csa_seed=prefix(:lenpre)//'.CSA.seed'
2059 csa_history=prefix(:lenpre)//'.CSA.history'
2060 csa_bank=prefix(:lenpre)//'.CSA.bank'
2061 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2062 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2063 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2064 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2065 csa_int=prefix(:lenpre)//'.int'
2066 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2067 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2068 csa_in=prefix(:lenpre)//'.CSA.in'
2069 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2072 write (iout,'(80(1h-))')
2073 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2074 write (iout,'(80(1h-))')
2075 write (iout,*) "Input file : ",
2076 & pref_orig(:ilen(pref_orig))//'.inp'
2077 write (iout,*) "Output file : ",
2078 & outname(:ilen(outname))
2080 write (iout,*) "Sidechain potential file : ",
2081 & sidename(:ilen(sidename))
2083 write (iout,*) "SCp potential file : ",
2084 & scpname(:ilen(scpname))
2086 write (iout,*) "Electrostatic potential file : ",
2087 & elename(:ilen(elename))
2088 write (iout,*) "Cumulant coefficient file : ",
2089 & fouriername(:ilen(fouriername))
2090 write (iout,*) "Torsional parameter file : ",
2091 & torname(:ilen(torname))
2092 write (iout,*) "Double torsional parameter file : ",
2093 & tordname(:ilen(tordname))
2094 write (iout,*) "SCCOR parameter file : ",
2095 & sccorname(:ilen(sccorname))
2096 write (iout,*) "Bond & inertia constant file : ",
2097 & bondname(:ilen(bondname))
2098 write (iout,*) "Bending parameter file : ",
2099 & thetname(:ilen(thetname))
2100 write (iout,*) "Rotamer parameter file : ",
2101 & rotname(:ilen(rotname))
2102 write (iout,*) "Thetpdb parameter file : ",
2103 & thetname_pdb(:ilen(thetname_pdb))
2104 write (iout,*) "Threading database : ",
2105 & patname(:ilen(patname))
2107 &write (iout,*)" DIRTMP : ",
2109 write (iout,'(80(1h-))')
2113 c----------------------------------------------------------------------------
2114 subroutine card_concat(card)
2115 implicit real*8 (a-h,o-z)
2116 include 'DIMENSIONS'
2117 include 'COMMON.IOUNITS'
2119 character*80 karta,ucase
2121 read (inp,'(a)') karta
2124 do while (karta(80:80).eq.'&')
2125 card=card(:ilen(card)+1)//karta(:79)
2126 read (inp,'(a)') karta
2129 card=card(:ilen(card)+1)//karta
2132 c----------------------------------------------------------------------------------
2134 implicit real*8 (a-h,o-z)
2135 include 'DIMENSIONS'
2136 include 'COMMON.CHAIN'
2137 include 'COMMON.IOUNITS'
2139 open(irest2,file=rest2name,status='unknown')
2140 read(irest2,*) totT,EK,potE,totE,t_bath
2142 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2145 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2148 read (irest2,*) iset
2153 c---------------------------------------------------------------------------------
2154 subroutine read_fragments
2155 implicit real*8 (a-h,o-z)
2156 include 'DIMENSIONS'
2160 include 'COMMON.SETUP'
2161 include 'COMMON.CHAIN'
2162 include 'COMMON.IOUNITS'
2164 include 'COMMON.CONTROL'
2165 read(inp,*) nset,nfrag,npair,nfrag_back
2166 if(me.eq.king.or..not.out1file)
2167 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2168 & " nfrag_back",nfrag_back
2170 read(inp,*) mset(iset)
2172 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2174 if(me.eq.king.or..not.out1file)
2175 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2176 & ifrag(2,i,iset), qinfrag(i,iset)
2179 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2181 if(me.eq.king.or..not.out1file)
2182 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2183 & ipair(2,i,iset), qinpair(i,iset)
2186 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2187 & wfrag_back(3,i,iset),
2188 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2189 if(me.eq.king.or..not.out1file)
2190 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2191 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2196 c-------------------------------------------------------------------------------
2197 subroutine read_dist_constr
2198 implicit real*8 (a-h,o-z)
2199 include 'DIMENSIONS'
2203 include 'COMMON.SETUP'
2204 include 'COMMON.CONTROL'
2205 include 'COMMON.CHAIN'
2206 include 'COMMON.IOUNITS'
2207 include 'COMMON.SBRIDGE'
2208 integer ifrag_(2,100),ipair_(2,100)
2209 double precision wfrag_(100),wpair_(100)
2210 character*500 controlcard
2211 c write (iout,*) "Calling read_dist_constr"
2212 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2214 call card_concat(controlcard)
2215 call readi(controlcard,"NFRAG",nfrag_,0)
2216 call readi(controlcard,"NPAIR",npair_,0)
2217 call readi(controlcard,"NDIST",ndist_,0)
2218 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2219 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2220 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2221 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2222 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2223 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2224 c write (iout,*) "IFRAG"
2226 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2228 c write (iout,*) "IPAIR"
2230 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2234 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2235 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2236 & ifrag_(2,i)=nstart_sup+nsup-1
2237 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2239 if (wfrag_(i).gt.0.0d0) then
2240 do j=ifrag_(1,i),ifrag_(2,i)-1
2241 do k=j+1,ifrag_(2,i)
2242 c write (iout,*) "j",j," k",k
2244 if (constr_dist.eq.1) then
2249 forcon(nhpb)=wfrag_(i)
2250 else if (constr_dist.eq.2) then
2251 if (ddjk.le.dist_cut) then
2256 forcon(nhpb)=wfrag_(i)
2263 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2266 if (.not.out1file .or. me.eq.king)
2267 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2268 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2270 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2271 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2278 if (wpair_(i).gt.0.0d0) then
2286 do j=ifrag_(1,ii),ifrag_(2,ii)
2287 do k=ifrag_(1,jj),ifrag_(2,jj)
2291 forcon(nhpb)=wpair_(i)
2292 dhpb(nhpb)=dist(j,k)
2294 if (.not.out1file .or. me.eq.king)
2295 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2296 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2298 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2299 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2306 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2307 if (forcon(nhpb+1).gt.0.0d0) then
2309 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2311 if (.not.out1file .or. me.eq.king)
2312 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2313 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2315 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2316 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2323 c-------------------------------------------------------------------------------
2325 subroutine flush(iu)
2330 subroutine flush(iu)
2335 c------------------------------------------------------------------------------
2336 subroutine copy_to_tmp(source)
2337 include "DIMENSIONS"
2338 include "COMMON.IOUNITS"
2339 character*(*) source
2340 character* 256 tmpfile
2344 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2345 inquire(file=tmpfile,exist=ex)
2347 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2348 & " to temporary directory..."
2349 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2350 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2354 c------------------------------------------------------------------------------
2355 subroutine move_from_tmp(source)
2356 include "DIMENSIONS"
2357 include "COMMON.IOUNITS"
2358 character*(*) source
2361 write (*,*) "Moving ",source(:ilen(source)),
2362 & " from temporary directory to working directory"
2363 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2364 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2367 c------------------------------------------------------------------------------
2368 subroutine random_init(seed)
2370 C Initialize random number generator
2372 implicit real*8 (a-h,o-z)
2373 include 'DIMENSIONS'
2376 logical OKRandom, prng_restart
2378 integer iseed_array(4)
2380 include 'COMMON.IOUNITS'
2381 include 'COMMON.TIME1'
2382 include 'COMMON.THREAD'
2383 include 'COMMON.SBRIDGE'
2384 include 'COMMON.CONTROL'
2385 include 'COMMON.MCM'
2386 include 'COMMON.MAP'
2387 include 'COMMON.HEADER'
2388 include 'COMMON.CSA'
2389 include 'COMMON.CHAIN'
2390 include 'COMMON.MUCA'
2392 include 'COMMON.FFIELD'
2393 include 'COMMON.SETUP'
2394 iseed=-dint(dabs(seed))
2395 if (iseed.eq.0) then
2396 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2397 & 'Random seed undefined. The program will stop.'
2398 write (*,'(/80(1h*)/20x,a/80(1h*))')
2399 & 'Random seed undefined. The program will stop.'
2401 call mpi_finalize(mpi_comm_world,ierr)
2403 stop 'Bad random seed.'
2406 if (fg_rank.eq.0) then
2410 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2411 OKRandom = prng_restart(me,iseed)
2414 tmp=65536.0d0**(4-i)
2415 iseed_array(i) = dint(seed/tmp)
2416 seed=seed-iseed_array(i)*tmp
2419 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2420 & (iseed_array(i),i=1,4)
2421 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2422 & (iseed_array(i),i=1,4)
2423 OKRandom = prng_restart(me,iseed_array)
2426 c r1 = prng_next(me)
2427 r1=ran_number(0.0D0,1.0D0)
2429 & write (iout,*) 'ran_num',r1
2430 if (r1.lt.0.0d0) OKRandom=.false.
2432 if (.not.OKRandom) then
2433 write (iout,*) 'PRNG IS NOT WORKING!!!'
2434 print *,'PRNG IS NOT WORKING!!!'
2437 call mpi_abort(mpi_comm_world,error_msg,ierr)
2440 write (iout,*) 'too many processors for parallel prng'
2441 write (*,*) 'too many processors for parallel prng'
2449 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)