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
348 mdpdb = index(controlcard,"MDPDB").gt.0
349 call reada(controlcard,"T_BATH",t_bath,300.0d0)
350 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
351 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
352 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
353 if (count_reset_moment.eq.0) count_reset_moment=1000000000
354 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
355 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
356 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
357 if (count_reset_vel.eq.0) count_reset_vel=1000000000
358 large = index(controlcard,"LARGE").gt.0
359 print_compon = index(controlcard,"PRINT_COMPON").gt.0
360 rattle = index(controlcard,"RATTLE").gt.0
361 c if performing umbrella sampling, fragments constrained are read from the fragment file
367 if(me.eq.king.or..not.out1file) then
369 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
371 write (iout,'(a)') "The units are:"
372 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
373 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
374 & " acceleration: angstrom/(48.9 fs)**2"
375 write (iout,'(a)') "energy: kcal/mol, temperature: K"
377 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
378 write (iout,'(a60,f10.5,a)')
379 & "Initial time step of numerical integration:",d_time,
381 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
383 write (iout,'(2a,i4,a)')
384 & "A-MTS algorithm used; initial time step for fast-varying",
385 & " short-range forces split into",ntime_split," steps."
386 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
387 & r_cut," lambda",rlamb
389 write (iout,'(2a,f10.5)')
390 & "Maximum acceleration threshold to reduce the time step",
391 & "/increase split number:",damax
392 write (iout,'(2a,f10.5)')
393 & "Maximum predicted energy drift to reduce the timestep",
394 & "/increase split number:",edriftmax
395 write (iout,'(a60,f10.5)')
396 & "Maximum velocity threshold to reduce velocities:",dvmax
397 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
398 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
399 if (rattle) write (iout,'(a60)')
400 & "Rattle algorithm used to constrain the virtual bonds"
404 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
405 call reada(controlcard,"RWAT",rwat,1.4d0)
406 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
407 surfarea=index(controlcard,"SURFAREA").gt.0
408 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
409 if(me.eq.king.or..not.out1file)then
410 write (iout,'(/a,$)') "Langevin dynamics calculation"
413 & " with direct integration of Langevin equations"
414 else if (lang.eq.2) then
415 write (iout,'(a/)') " with TINKER stochasic MD integrator"
416 else if (lang.eq.3) then
417 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
418 else if (lang.eq.4) then
419 write (iout,'(a/)') " in overdamped mode"
421 write (iout,'(//a,i5)')
422 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
425 write (iout,'(a60,f10.5)') "Temperature:",t_bath
426 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
427 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
428 write (iout,'(a60,f10.5)')
429 & "Scaling factor of the friction forces:",scal_fric
430 if (surfarea) write (iout,'(2a,i10,a)')
431 & "Friction coefficients will be scaled by solvent-accessible",
432 & " surface area every",reset_fricmat," steps."
434 c Calculate friction coefficients and bounds of stochastic forces
435 eta=6*pi*cPoise*etawat
436 if(me.eq.king.or..not.out1file)
437 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
439 gamp=scal_fric*(pstok+rwat)*eta
440 stdfp=dsqrt(2*Rb*t_bath/d_time)
442 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
443 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
445 if(me.eq.king.or..not.out1file)then
446 write (iout,'(/2a/)')
447 & "Radii of site types and friction coefficients and std's of",
448 & " stochastic forces of fully exposed sites"
449 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
451 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
452 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
456 if(me.eq.king.or..not.out1file)then
457 write (iout,'(a)') "Berendsen bath calculation"
458 write (iout,'(a60,f10.5)') "Temperature:",t_bath
459 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
461 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
462 & count_reset_moment," steps"
464 & write (iout,'(a,i10,a)')
465 & "Velocities will be reset at random every",count_reset_vel,
469 if(me.eq.king.or..not.out1file)
470 & write (iout,'(a31)') "Microcanonical mode calculation"
472 if(me.eq.king.or..not.out1file)then
473 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
475 write(iout,*) "MD running with constraints."
476 write(iout,*) "Equilibration time ", eq_time, " mtus."
477 write(iout,*) "Constraining ", nfrag," fragments."
478 write(iout,*) "Length of each fragment, weight and q0:"
480 write (iout,*) "Set of restraints #",iset
482 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
483 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
485 write(iout,*) "constraints between ", npair, "fragments."
486 write(iout,*) "constraint pairs, weights and q0:"
488 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
489 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
491 write(iout,*) "angle constraints within ", nfrag_back,
492 & "backbone fragments."
493 write(iout,*) "fragment, weights:"
495 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
496 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
497 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
500 iset=mod(kolor,nset)+1
503 if(me.eq.king.or..not.out1file)
504 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
507 c------------------------------------------------------------------------------
510 C Read molecular data.
512 implicit real*8 (a-h,o-z)
518 include 'COMMON.IOUNITS'
521 include 'COMMON.INTERACT'
522 include 'COMMON.LOCAL'
523 include 'COMMON.NAMES'
524 include 'COMMON.CHAIN'
525 include 'COMMON.FFIELD'
526 include 'COMMON.SBRIDGE'
527 include 'COMMON.HEADER'
528 include 'COMMON.CONTROL'
529 include 'COMMON.DBASE'
530 include 'COMMON.THREAD'
531 include 'COMMON.CONTACTS'
532 include 'COMMON.TORCNSTR'
533 include 'COMMON.TIME1'
534 include 'COMMON.BOUNDS'
536 include 'COMMON.SETUP'
537 character*4 sequence(maxres)
539 double precision x(maxvar)
540 character*256 pdbfile
541 character*320 weightcard
542 character*80 weightcard_t,ucase
543 dimension itype_pdb(maxres)
544 common /pizda/ itype_pdb
545 logical seq_comp,fail
546 double precision energia(0:n_ene)
552 C Read weights of the subsequent energy terms.
553 call card_concat(weightcard)
554 call reada(weightcard,'WLONG',wlong,1.0D0)
555 call reada(weightcard,'WSC',wsc,wlong)
556 call reada(weightcard,'WSCP',wscp,wlong)
557 call reada(weightcard,'WELEC',welec,1.0D0)
558 call reada(weightcard,'WVDWPP',wvdwpp,welec)
559 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
560 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
561 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
562 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
563 call reada(weightcard,'WTURN3',wturn3,1.0D0)
564 call reada(weightcard,'WTURN4',wturn4,1.0D0)
565 call reada(weightcard,'WTURN6',wturn6,1.0D0)
566 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
567 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
568 call reada(weightcard,'WBOND',wbond,1.0D0)
569 call reada(weightcard,'WTOR',wtor,1.0D0)
570 call reada(weightcard,'WTORD',wtor_d,1.0D0)
571 call reada(weightcard,'WANG',wang,1.0D0)
572 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
573 call reada(weightcard,'SCAL14',scal14,0.4D0)
574 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
575 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
576 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
577 call reada(weightcard,'TEMP0',temp0,300.0d0)
578 if (index(weightcard,'SOFT').gt.0) ipot=6
579 C 12/1/95 Added weight for the multi-body term WCORR
580 call reada(weightcard,'WCORRH',wcorr,1.0D0)
581 if (wcorr4.gt.0.0d0) wcorr=wcorr4
601 if(me.eq.king.or..not.out1file)
602 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
603 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
605 10 format (/'Energy-term weights (unscaled):'//
606 & 'WSCC= ',f10.6,' (SC-SC)'/
607 & 'WSCP= ',f10.6,' (SC-p)'/
608 & 'WELEC= ',f10.6,' (p-p electr)'/
609 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
610 & 'WBOND= ',f10.6,' (stretching)'/
611 & 'WANG= ',f10.6,' (bending)'/
612 & 'WSCLOC= ',f10.6,' (SC local)'/
613 & 'WTOR= ',f10.6,' (torsional)'/
614 & 'WTORD= ',f10.6,' (double torsional)'/
615 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
616 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
617 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
618 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
619 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
620 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
621 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
622 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
623 & 'WTURN6= ',f10.6,' (turns, 6th order)')
624 if(me.eq.king.or..not.out1file)then
625 if (wcorr4.gt.0.0d0) then
626 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
627 & 'between contact pairs of peptide groups'
628 write (iout,'(2(a,f5.3/))')
629 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
630 & 'Range of quenching the correlation terms:',2*delt_corr
631 else if (wcorr.gt.0.0d0) then
632 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
633 & 'between contact pairs of peptide groups'
635 write (iout,'(a,f8.3)')
636 & 'Scaling factor of 1,4 SC-p interactions:',scal14
637 write (iout,'(a,f8.3)')
638 & 'General scaling factor of SC-p interactions:',scalscp
640 r0_corr=cutoff_corr-delt_corr
642 aad(i,1)=scalscp*aad(i,1)
643 aad(i,2)=scalscp*aad(i,2)
644 bad(i,1)=scalscp*bad(i,1)
645 bad(i,2)=scalscp*bad(i,2)
647 call rescale_weights(t_bath)
648 if(me.eq.king.or..not.out1file)
649 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
650 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
652 22 format (/'Energy-term weights (scaled):'//
653 & 'WSCC= ',f10.6,' (SC-SC)'/
654 & 'WSCP= ',f10.6,' (SC-p)'/
655 & 'WELEC= ',f10.6,' (p-p electr)'/
656 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
657 & 'WBOND= ',f10.6,' (stretching)'/
658 & 'WANG= ',f10.6,' (bending)'/
659 & 'WSCLOC= ',f10.6,' (SC local)'/
660 & 'WTOR= ',f10.6,' (torsional)'/
661 & 'WTORD= ',f10.6,' (double torsional)'/
662 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
663 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
664 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
665 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
666 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
667 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
668 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
669 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
670 & 'WTURN6= ',f10.6,' (turns, 6th order)')
671 if(me.eq.king.or..not.out1file)
672 & write (iout,*) "Reference temperature for weights calculation:",
674 call reada(weightcard,"D0CM",d0cm,3.78d0)
675 call reada(weightcard,"AKCM",akcm,15.1d0)
676 call reada(weightcard,"AKTH",akth,11.0d0)
677 call reada(weightcard,"AKCT",akct,12.0d0)
678 call reada(weightcard,"V1SS",v1ss,-1.08d0)
679 call reada(weightcard,"V2SS",v2ss,7.61d0)
680 call reada(weightcard,"V3SS",v3ss,13.7d0)
681 call reada(weightcard,"EBR",ebr,-5.50D0)
682 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
684 dyn_ss_mask(i)=.false.
688 dyn_ssbond_ij(i,j)=1.0d300
691 call reada(weightcard,"HT",Ht,0.0D0)
693 ss_depth=ebr/wsc-0.25*eps(1,1)
694 Ht=Ht/wsc-0.25*eps(1,1)
695 akcm=akcm*wstrain/wsc
696 akth=akth*wstrain/wsc
697 akct=akct*wstrain/wsc
698 v1ss=v1ss*wstrain/wsc
699 v2ss=v2ss*wstrain/wsc
700 v3ss=v3ss*wstrain/wsc
702 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
705 if(me.eq.king.or..not.out1file) then
706 write (iout,*) "Parameters of the SS-bond potential:"
707 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
709 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
710 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
711 write (iout,*)" HT",Ht
712 print *,'indpdb=',indpdb,' pdbref=',pdbref
714 if (indpdb.gt.0 .or. pdbref) then
715 read(inp,'(a)') pdbfile
716 if(me.eq.king.or..not.out1file)
717 & write (iout,'(2a)') 'PDB data will be read from file ',
718 & pdbfile(:ilen(pdbfile))
719 open(ipdbin,file=pdbfile,status='old',err=33)
721 33 write (iout,'(a)') 'Error opening PDB file.'
724 c print *,'Begin reading pdb data'
726 c print *,'Finished reading pdb data'
727 if(me.eq.king.or..not.out1file)
728 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
729 & ' nstart_sup=',nstart_sup
731 itype_pdb(i)=itype(i)
735 nct=nstart_sup+nsup-1
736 call contact(.false.,ncont_ref,icont_ref,co)
739 if(me.eq.king.or..not.out1file)
740 & write(iout,*)'Adding sidechains'
744 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
747 do while (fail.and.nsi.le.maxsi)
748 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
751 if(fail) write(iout,*)'Adding sidechain failed for res ',
752 & i,' after ',nsi,' trials'
757 if (indpdb.eq.0) then
758 C Read sequence if not taken from the pdb file.
760 c print *,'nres=',nres
761 if (iscode.gt.0) then
762 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
764 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
766 C Convert sequence to numeric code
768 itype(i)=rescode(i,sequence(i),iscode)
770 C Assign initial virtual bond lengths
776 vbld(i+nres)=dsc(iabs(itype(i)))
777 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
778 c write (iout,*) "i",i," itype",itype(i),
779 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
783 c print '(20i4)',(itype(i),i=1,nres)
786 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
788 if (itype(i).eq.ntyp1) then
792 else if (iabs(itype(i+1)).ne.20) then
794 else if (iabs(itype(i)).ne.20) then
801 if(me.eq.king.or..not.out1file)then
802 write (iout,*) "ITEL"
804 write (iout,*) i,itype(i),itel(i)
806 print *,'Call Read_Bridge.'
809 C 8/13/98 Set limits to generating the dihedral angles
814 read (inp,*) ndih_constr
815 if (ndih_constr.gt.0) then
817 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
818 if(me.eq.king.or..not.out1file)then
820 & 'There are',ndih_constr,' constraints on phi angles.'
822 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
826 phi0(i)=deg2rad*phi0(i)
827 drange(i)=deg2rad*drange(i)
829 if(me.eq.king.or..not.out1file)
830 & write (iout,*) 'FTORS',ftors
833 phibound(1,ii) = phi0(i)-drange(i)
834 phibound(2,ii) = phi0(i)+drange(i)
841 write (iout,'(a)') 'Boundaries in phi angle sampling:'
843 write (iout,'(a3,i5,2f10.1)')
844 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
850 cd print *,'NNT=',NNT,' NCT=',NCT
851 if (itype(1).eq.ntyp1) nnt=2
852 if (itype(nres).eq.ntyp1) nct=nct-1
854 if(me.eq.king.or..not.out1file)
855 & write (iout,'(a,i3)') 'nsup=',nsup
857 if (nsup.le.(nct-nnt+1)) then
858 do i=0,nct-nnt+1-nsup
859 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
865 & 'Error - sequences to be superposed do not match.'
868 do i=0,nsup-(nct-nnt+1)
869 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
871 nstart_sup=nstart_sup+i
877 & 'Error - sequences to be superposed do not match.'
880 if (nsup.eq.0) nsup=nct-nnt
881 if (nstart_sup.eq.0) nstart_sup=nnt
882 if (nstart_seq.eq.0) nstart_seq=nnt
883 if(me.eq.king.or..not.out1file)
884 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
885 & ' nstart_seq=',nstart_seq
887 c--- Zscore rms -------
888 if (nz_start.eq.0) nz_start=nnt
889 if (nz_end.eq.0 .and. nsup.gt.0) then
891 else if (nz_end.eq.0) then
894 if(me.eq.king.or..not.out1file)then
895 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
896 write (iout,*) 'IZ_SC=',iz_sc
898 c----------------------
901 if (.not.pdbref) then
902 call read_angles(inp,*38)
904 38 write (iout,'(a)') 'Error reading reference structure.'
906 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
907 stop 'Error reading reference structure'
911 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
921 call contact(.true.,ncont_ref,icont_ref,co)
923 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
925 if (constr_dist.gt.0) call read_dist_constr
926 write (iout,*) "After read_dist_constr nhpb",nhpb
928 if(me.eq.king.or..not.out1file)
929 & write (iout,*) 'Contact order:',co
931 if(me.eq.king.or..not.out1file)
932 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
935 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
937 if(me.eq.king.or..not.out1file)
938 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
939 & icont_ref(1,i),' ',
940 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
944 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
945 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
946 & modecalc.ne.10) then
947 C If input structure hasn't been supplied from the PDB file read or generate
949 if (iranconf.eq.0 .and. .not. extconf) then
950 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
951 & write (iout,'(a)') 'Initial geometry will be read in.'
953 read(inp,'(8f10.5)',end=36,err=36)
954 & ((c(l,k),l=1,3),k=1,nres),
955 & ((c(l,k+nres),l=1,3),k=nnt,nct)
956 write (iout,*) "Exit READ_CART"
957 write (iout,'(8f10.5)')
958 & ((c(l,k),l=1,3),k=1,nres),
959 & ((c(l,k+nres),l=1,3),k=nnt,nct)
960 call int_from_cart1(.true.)
961 write (iout,*) "Finish INT_TO_CART"
964 dc(j,i)=c(j,i+1)-c(j,i)
965 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
969 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
971 dc(j,i+nres)=c(j,i+nres)-c(j,i)
972 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
978 call read_angles(inp,*36)
981 36 write (iout,'(a)') 'Error reading angle file.'
983 call mpi_finalize( MPI_COMM_WORLD,IERR )
985 stop 'Error reading angle file.'
987 else if (extconf) then
988 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
989 & write (iout,'(a)') 'Extended chain initial geometry.'
991 theta(i)=90d0*deg2rad
997 alph(i)=110d0*deg2rad
1000 omeg(i)=-120d0*deg2rad
1001 if (itype(i).le.0) omeg(i)=-omeg(i)
1004 if(me.eq.king.or..not.out1file)
1005 & write (iout,'(a)') 'Random-generated initial geometry.'
1009 if (me.eq.king .or. fg_rank.eq.0 .and. (
1010 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1014 call gen_rand_conf(itmp,*30)
1016 30 write (iout,*) 'Failed to generate random conformation',
1017 & ', itrial=',itrial
1018 write (*,*) 'Processor:',me,
1019 & ' Failed to generate random conformation',
1029 write (iout,'(a,i3,a)') 'Processor:',me,
1030 & ' error in generating random conformation.'
1031 write (*,'(a,i3,a)') 'Processor:',me,
1032 & ' error in generating random conformation.'
1035 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1041 call gen_rand_conf(itmp,*30)
1043 30 write (iout,*) 'Failed to generate random conformation',
1044 & ', itrial=',itrial
1045 write (*,*) 'Failed to generate random conformation',
1046 & ', itrial=',itrial
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.'
1056 elseif (modecalc.eq.4) then
1057 read (inp,'(a)') intinname
1058 open (intin,file=intinname,status='old',err=333)
1059 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1060 & write (iout,'(a)') 'intinname',intinname
1061 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1063 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1065 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1067 stop 'Error opening angle file.'
1071 C Generate distance constraints, if the PDB structure is to be regularized.
1072 if (nthread.gt.0) then
1073 call read_threadbase
1076 if (me.eq.king .or. .not. out1file)
1078 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1079 write (iout,'(/a,i3,a)')
1080 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1081 write (iout,'(20i4)') (iss(i),i=1,ns)
1083 write(iout,*)"Running with dynamic disulfide-bond formation"
1085 write (iout,'(/a/)') 'Pre-formed links are:'
1091 if (me.eq.king.or..not.out1file)
1092 & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1093 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1099 if (ns.gt.0.and.dyn_ss) then
1103 forcon(i-nss)=forcon(i)
1110 dyn_ss_mask(iss(i))=.true.
1113 if (i2ndstr.gt.0) call secstrp2dihc
1114 c call geom_to_var(nvar,x)
1115 c call etotal(energia(0))
1116 c call enerprint(energia(0))
1117 c call briefout(0,etot)
1119 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1120 cd write (iout,'(a)') 'Variable list:'
1121 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1123 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1124 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1125 & 'Processor',myrank,': end reading molecular data.'
1129 c--------------------------------------------------------------------------
1130 logical function seq_comp(itypea,itypeb,length)
1132 integer length,itypea(length),itypeb(length)
1135 if (itypea(i).ne.itypeb(i)) then
1143 c-----------------------------------------------------------------------------
1144 subroutine read_bridge
1145 C Read information about disulfide bridges.
1146 implicit real*8 (a-h,o-z)
1147 include 'DIMENSIONS'
1151 include 'COMMON.IOUNITS'
1152 include 'COMMON.GEO'
1153 include 'COMMON.VAR'
1154 include 'COMMON.INTERACT'
1155 include 'COMMON.LOCAL'
1156 include 'COMMON.NAMES'
1157 include 'COMMON.CHAIN'
1158 include 'COMMON.FFIELD'
1159 include 'COMMON.SBRIDGE'
1160 include 'COMMON.HEADER'
1161 include 'COMMON.CONTROL'
1162 include 'COMMON.DBASE'
1163 include 'COMMON.THREAD'
1164 include 'COMMON.TIME1'
1165 include 'COMMON.SETUP'
1166 C Read bridging residues.
1167 read (inp,*) ns,(iss(i),i=1,ns)
1169 if(me.eq.king.or..not.out1file)
1170 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1171 C Check whether the specified bridging residues are cystines.
1173 if (itype(iss(i)).ne.1) then
1174 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1175 & 'Do you REALLY think that the residue ',
1176 & restyp(itype(iss(i))),i,
1177 & ' can form a disulfide bridge?!!!'
1178 write (*,'(2a,i3,a)')
1179 & 'Do you REALLY think that the residue ',
1180 & restyp(itype(iss(i))),i,
1181 & ' can form a disulfide bridge?!!!'
1183 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1188 C Read preformed bridges.
1190 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1192 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1195 C Check if the residues involved in bridges are in the specified list of
1196 C bridging residues.
1199 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1200 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1201 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1202 & ' contains residues present in other pairs.'
1203 write (*,'(a,i3,a)') 'Disulfide pair',i,
1204 & ' contains residues present in other pairs.'
1206 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1212 if (ihpb(i).eq.iss(j)) goto 10
1214 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1217 if (jhpb(i).eq.iss(j)) goto 20
1219 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1225 ihpb(i)=ihpb(i)+nres
1226 jhpb(i)=jhpb(i)+nres
1232 c----------------------------------------------------------------------------
1233 subroutine read_x(kanal,*)
1234 implicit real*8 (a-h,o-z)
1235 include 'DIMENSIONS'
1236 include 'COMMON.GEO'
1237 include 'COMMON.VAR'
1238 include 'COMMON.CHAIN'
1239 include 'COMMON.IOUNITS'
1240 include 'COMMON.CONTROL'
1241 include 'COMMON.LOCAL'
1242 include 'COMMON.INTERACT'
1243 c Read coordinates from input
1245 read(kanal,'(8f10.5)',end=10,err=10)
1246 & ((c(l,k),l=1,3),k=1,nres),
1247 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1250 c(j,2*nres)=c(j,nres)
1252 call int_from_cart1(.false.)
1255 dc(j,i)=c(j,i+1)-c(j,i)
1256 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1260 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1262 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1263 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1271 c----------------------------------------------------------------------------
1272 subroutine read_threadbase
1273 implicit real*8 (a-h,o-z)
1274 include 'DIMENSIONS'
1275 include 'COMMON.IOUNITS'
1276 include 'COMMON.GEO'
1277 include 'COMMON.VAR'
1278 include 'COMMON.INTERACT'
1279 include 'COMMON.LOCAL'
1280 include 'COMMON.NAMES'
1281 include 'COMMON.CHAIN'
1282 include 'COMMON.FFIELD'
1283 include 'COMMON.SBRIDGE'
1284 include 'COMMON.HEADER'
1285 include 'COMMON.CONTROL'
1286 include 'COMMON.DBASE'
1287 include 'COMMON.THREAD'
1288 include 'COMMON.TIME1'
1289 C Read pattern database for threading.
1290 read (icbase,*) nseq
1292 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1293 & nres_base(2,i),nres_base(3,i)
1294 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1296 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1297 c & nres_base(2,i),nres_base(3,i)
1298 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1302 if (weidis.eq.0.0D0) weidis=0.1D0
1311 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1312 write (iout,'(a,i5)') 'nexcl: ',nexcl
1313 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1316 c------------------------------------------------------------------------------
1317 subroutine setup_var
1318 implicit real*8 (a-h,o-z)
1319 include 'DIMENSIONS'
1320 include 'COMMON.IOUNITS'
1321 include 'COMMON.GEO'
1322 include 'COMMON.VAR'
1323 include 'COMMON.INTERACT'
1324 include 'COMMON.LOCAL'
1325 include 'COMMON.NAMES'
1326 include 'COMMON.CHAIN'
1327 include 'COMMON.FFIELD'
1328 include 'COMMON.SBRIDGE'
1329 include 'COMMON.HEADER'
1330 include 'COMMON.CONTROL'
1331 include 'COMMON.DBASE'
1332 include 'COMMON.THREAD'
1333 include 'COMMON.TIME1'
1334 C Set up variable list.
1340 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1342 ialph(i,1)=nvar+nside
1346 if (indphi.gt.0) then
1348 else if (indback.gt.0) then
1353 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1356 c----------------------------------------------------------------------------
1357 subroutine gen_dist_constr
1358 C Generate CA distance constraints.
1359 implicit real*8 (a-h,o-z)
1360 include 'DIMENSIONS'
1361 include 'COMMON.IOUNITS'
1362 include 'COMMON.GEO'
1363 include 'COMMON.VAR'
1364 include 'COMMON.INTERACT'
1365 include 'COMMON.LOCAL'
1366 include 'COMMON.NAMES'
1367 include 'COMMON.CHAIN'
1368 include 'COMMON.FFIELD'
1369 include 'COMMON.SBRIDGE'
1370 include 'COMMON.HEADER'
1371 include 'COMMON.CONTROL'
1372 include 'COMMON.DBASE'
1373 include 'COMMON.THREAD'
1374 include 'COMMON.TIME1'
1375 dimension itype_pdb(maxres)
1376 common /pizda/ itype_pdb
1378 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1379 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1380 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1382 do i=nstart_sup,nstart_sup+nsup-1
1383 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1384 cd & ' seq_pdb', restyp(itype_pdb(i))
1385 do j=i+2,nstart_sup+nsup-1
1387 ihpb(nhpb)=i+nstart_seq-nstart_sup
1388 jhpb(nhpb)=j+nstart_seq-nstart_sup
1390 dhpb(nhpb)=dist(i,j)
1393 cd write (iout,'(a)') 'Distance constraints:'
1398 cd if (ii.gt.nres) then
1403 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1404 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1405 cd & dhpb(i),forcon(i)
1409 c----------------------------------------------------------------------------
1411 implicit real*8 (a-h,o-z)
1412 include 'DIMENSIONS'
1413 include 'COMMON.MAP'
1414 include 'COMMON.IOUNITS'
1415 character*3 angid(4) /'THE','PHI','ALP','OME'/
1416 character*80 mapcard,ucase
1418 read (inp,'(a)') mapcard
1419 mapcard=ucase(mapcard)
1420 if (index(mapcard,'PHI').gt.0) then
1422 else if (index(mapcard,'THE').gt.0) then
1424 else if (index(mapcard,'ALP').gt.0) then
1426 else if (index(mapcard,'OME').gt.0) then
1429 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1430 stop 'Error - illegal variable spec in MAP card.'
1432 call readi (mapcard,'RES1',res1(imap),0)
1433 call readi (mapcard,'RES2',res2(imap),0)
1434 if (res1(imap).eq.0) then
1435 res1(imap)=res2(imap)
1436 else if (res2(imap).eq.0) then
1437 res2(imap)=res1(imap)
1439 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1441 & 'Error - illegal definition of variable group in MAP.'
1442 stop 'Error - illegal definition of variable group in MAP.'
1444 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1445 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1446 call readi(mapcard,'NSTEP',nstep(imap),0)
1447 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1449 & 'Illegal boundary and/or step size specification in MAP.'
1450 stop 'Illegal boundary and/or step size specification in MAP.'
1455 c----------------------------------------------------------------------------
1457 implicit real*8 (a-h,o-z)
1458 include 'DIMENSIONS'
1459 include 'COMMON.IOUNITS'
1460 include 'COMMON.GEO'
1461 include 'COMMON.CSA'
1462 include 'COMMON.BANK'
1463 include 'COMMON.CONTROL'
1465 character*620 mcmcard
1466 call card_concat(mcmcard)
1468 call readi(mcmcard,'NCONF',nconf,50)
1469 call readi(mcmcard,'NADD',nadd,0)
1470 call readi(mcmcard,'JSTART',jstart,1)
1471 call readi(mcmcard,'JEND',jend,1)
1472 call readi(mcmcard,'NSTMAX',nstmax,500000)
1473 call readi(mcmcard,'N0',n0,1)
1474 call readi(mcmcard,'N1',n1,6)
1475 call readi(mcmcard,'N2',n2,4)
1476 call readi(mcmcard,'N3',n3,0)
1477 call readi(mcmcard,'N4',n4,0)
1478 call readi(mcmcard,'N5',n5,0)
1479 call readi(mcmcard,'N6',n6,10)
1480 call readi(mcmcard,'N7',n7,0)
1481 call readi(mcmcard,'N8',n8,0)
1482 call readi(mcmcard,'N9',n9,0)
1483 call readi(mcmcard,'N14',n14,0)
1484 call readi(mcmcard,'N15',n15,0)
1485 call readi(mcmcard,'N16',n16,0)
1486 call readi(mcmcard,'N17',n17,0)
1487 call readi(mcmcard,'N18',n18,0)
1489 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1491 call readi(mcmcard,'NDIFF',ndiff,2)
1492 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1493 call readi(mcmcard,'IS1',is1,1)
1494 call readi(mcmcard,'IS2',is2,8)
1495 call readi(mcmcard,'NRAN0',nran0,4)
1496 call readi(mcmcard,'NRAN1',nran1,2)
1497 call readi(mcmcard,'IRR',irr,1)
1498 call readi(mcmcard,'NSEED',nseed,20)
1499 call readi(mcmcard,'NTOTAL',ntotal,10000)
1500 call reada(mcmcard,'CUT1',cut1,2.0d0)
1501 call reada(mcmcard,'CUT2',cut2,5.0d0)
1502 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1503 call readi(mcmcard,'ICMAX',icmax,3)
1504 call readi(mcmcard,'IRESTART',irestart,0)
1505 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1508 call reada(mcmcard,'DELE',dele,20.0d0)
1509 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1510 call readi(mcmcard,'IREF',iref,0)
1511 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1512 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1513 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1514 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1515 write (iout,*) "NCONF_IN",nconf_in
1518 c----------------------------------------------------------------------------
1519 cfmc subroutine mcmfread
1520 cfmc implicit real*8 (a-h,o-z)
1521 cfmc include 'DIMENSIONS'
1522 cfmc include 'COMMON.MCMF'
1523 cfmc include 'COMMON.IOUNITS'
1524 cfmc include 'COMMON.GEO'
1525 cfmc character*80 ucase
1526 cfmc character*620 mcmcard
1527 cfmc call card_concat(mcmcard)
1529 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1530 cfmc write(iout,*)'MAXRANT=',maxrant
1531 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1532 cfmc write(iout,*)'MAXFAM=',maxfam
1533 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1534 cfmc write(iout,*)'NNET1=',nnet1
1535 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1536 cfmc write(iout,*)'NNET2=',nnet2
1537 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1538 cfmc write(iout,*)'NNET3=',nnet3
1539 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1540 cfmc write(iout,*)'ILASTT=',ilastt
1541 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1542 cfmc write(iout,*)'MAXSTR=',maxstr
1543 cfmc maxstr_f=maxstr/maxfam
1544 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1545 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1546 cfmc write(iout,*)'NMCMF=',nmcmf
1547 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1548 cfmc write(iout,*)'IFOCUS=',ifocus
1549 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1550 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1551 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1552 cfmc write(iout,*)'INTPRT=',intprt
1553 cfmc call readi(mcmcard,'IPRT',iprt,100)
1554 cfmc write(iout,*)'IPRT=',iprt
1555 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1556 cfmc write(iout,*)'IMAXTR=',imaxtr
1557 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1558 cfmc write(iout,*)'MAXEVEN=',maxeven
1559 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1560 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1561 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1562 cfmc write(iout,*)'INIMIN=',inimin
1563 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1564 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1565 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1566 cfmc write(iout,*)'NTHREAD=',nthread
1567 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1568 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1569 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1570 cfmc write(iout,*)'MAXPERT=',maxpert
1571 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1572 cfmc write(iout,*)'IRMSD=',irmsd
1573 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1574 cfmc write(iout,*)'DENEMIN=',denemin
1575 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1576 cfmc write(iout,*)'RCUT1S=',rcut1s
1577 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1578 cfmc write(iout,*)'RCUT1E=',rcut1e
1579 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1580 cfmc write(iout,*)'RCUT2S=',rcut2s
1581 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1582 cfmc write(iout,*)'RCUT2E=',rcut2e
1583 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1584 cfmc write(iout,*)'DPERT1=',d_pert1
1585 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1586 cfmc write(iout,*)'DPERT1A=',d_pert1a
1587 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1588 cfmc write(iout,*)'DPERT2=',d_pert2
1589 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1590 cfmc write(iout,*)'DPERT2A=',d_pert2a
1591 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1592 cfmc write(iout,*)'DPERT2B=',d_pert2b
1593 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1594 cfmc write(iout,*)'DPERT2C=',d_pert2c
1595 cfmc d_pert1=deg2rad*d_pert1
1596 cfmc d_pert1a=deg2rad*d_pert1a
1597 cfmc d_pert2=deg2rad*d_pert2
1598 cfmc d_pert2a=deg2rad*d_pert2a
1599 cfmc d_pert2b=deg2rad*d_pert2b
1600 cfmc d_pert2c=deg2rad*d_pert2c
1601 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1602 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1603 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1604 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1605 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1606 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1607 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1608 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1609 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1610 cfmc write(iout,*)'RCUTINI=',rcutini
1611 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1612 cfmc write(iout,*)'GRAT=',grat
1613 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1614 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1618 c----------------------------------------------------------------------------
1620 implicit real*8 (a-h,o-z)
1621 include 'DIMENSIONS'
1622 include 'COMMON.MCM'
1623 include 'COMMON.MCE'
1624 include 'COMMON.IOUNITS'
1626 character*320 mcmcard
1627 call card_concat(mcmcard)
1628 call readi(mcmcard,'MAXACC',maxacc,100)
1629 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1630 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1631 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1632 call readi(mcmcard,'MAXREPM',maxrepm,200)
1633 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1634 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1635 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1636 call reada(mcmcard,'E_UP',e_up,5.0D0)
1637 call reada(mcmcard,'DELTE',delte,0.1D0)
1638 call readi(mcmcard,'NSWEEP',nsweep,5)
1639 call readi(mcmcard,'NSTEPH',nsteph,0)
1640 call readi(mcmcard,'NSTEPC',nstepc,0)
1641 call reada(mcmcard,'TMIN',tmin,298.0D0)
1642 call reada(mcmcard,'TMAX',tmax,298.0D0)
1643 call readi(mcmcard,'NWINDOW',nwindow,0)
1644 call readi(mcmcard,'PRINT_MC',print_mc,0)
1645 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1646 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1647 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1648 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1649 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1650 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1651 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1652 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1653 if (nwindow.gt.0) then
1654 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1656 winlen(i)=winend(i)-winstart(i)+1
1659 if (tmax.lt.tmin) tmax=tmin
1660 if (tmax.eq.tmin) then
1664 if (nstepc.gt.0 .and. nsteph.gt.0) then
1665 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1666 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1668 C Probabilities of different move types
1669 sumpro_type(0)=0.0D0
1670 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1671 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1672 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1673 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1674 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1675 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1676 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1678 print *,'i',i,' sumprotype',sumpro_type(i)
1679 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1680 print *,'i',i,' sumprotype',sumpro_type(i)
1684 c----------------------------------------------------------------------------
1685 subroutine read_minim
1686 implicit real*8 (a-h,o-z)
1687 include 'DIMENSIONS'
1688 include 'COMMON.MINIM'
1689 include 'COMMON.IOUNITS'
1691 character*320 minimcard
1692 call card_concat(minimcard)
1693 call readi(minimcard,'MAXMIN',maxmin,2000)
1694 call readi(minimcard,'MAXFUN',maxfun,5000)
1695 call readi(minimcard,'MINMIN',minmin,maxmin)
1696 call readi(minimcard,'MINFUN',minfun,maxmin)
1697 call reada(minimcard,'TOLF',tolf,1.0D-2)
1698 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1699 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1700 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1701 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1702 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1703 & 'Options in energy minimization:'
1704 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1705 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1706 & 'MinMin:',MinMin,' MinFun:',MinFun,
1707 & ' TolF:',TolF,' RTolF:',RTolF
1710 c----------------------------------------------------------------------------
1711 subroutine read_angles(kanal,*)
1712 implicit real*8 (a-h,o-z)
1713 include 'DIMENSIONS'
1714 include 'COMMON.GEO'
1715 include 'COMMON.VAR'
1716 include 'COMMON.CHAIN'
1717 include 'COMMON.IOUNITS'
1718 include 'COMMON.CONTROL'
1719 c Read angles from input
1721 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1722 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1723 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1724 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1727 c 9/7/01 avoid 180 deg valence angle
1728 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1730 theta(i)=deg2rad*theta(i)
1731 phi(i)=deg2rad*phi(i)
1732 alph(i)=deg2rad*alph(i)
1733 omeg(i)=deg2rad*omeg(i)
1738 c----------------------------------------------------------------------------
1739 subroutine reada(rekord,lancuch,wartosc,default)
1741 character*(*) rekord,lancuch
1742 double precision wartosc,default
1745 iread=index(rekord,lancuch)
1746 if (iread.eq.0) then
1750 iread=iread+ilen(lancuch)+1
1751 read (rekord(iread:),*,err=10,end=10) wartosc
1756 c----------------------------------------------------------------------------
1757 subroutine readi(rekord,lancuch,wartosc,default)
1759 character*(*) rekord,lancuch
1760 integer wartosc,default
1763 iread=index(rekord,lancuch)
1764 if (iread.eq.0) then
1768 iread=iread+ilen(lancuch)+1
1769 read (rekord(iread:),*,err=10,end=10) wartosc
1774 c----------------------------------------------------------------------------
1775 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1778 integer tablica(dim),default
1779 character*(*) rekord,lancuch
1786 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1787 if (iread.eq.0) return
1788 iread=iread+ilen(lancuch)+1
1789 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1792 c----------------------------------------------------------------------------
1793 subroutine multreada(rekord,lancuch,tablica,dim,default)
1796 double precision tablica(dim),default
1797 character*(*) rekord,lancuch
1804 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1805 if (iread.eq.0) return
1806 iread=iread+ilen(lancuch)+1
1807 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1810 c----------------------------------------------------------------------------
1811 subroutine openunits
1812 implicit real*8 (a-h,o-z)
1813 include 'DIMENSIONS'
1816 character*16 form,nodename
1819 include 'COMMON.SETUP'
1820 include 'COMMON.IOUNITS'
1822 include 'COMMON.CONTROL'
1823 integer lenpre,lenpot,ilen,lentmp
1825 character*3 out1file_text,ucase
1828 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1829 call getenv_loc("PREFIX",prefix)
1831 call getenv_loc("POT",pot)
1832 call getenv_loc("DIRTMP",tmpdir)
1833 call getenv_loc("CURDIR",curdir)
1834 call getenv_loc("OUT1FILE",out1file_text)
1835 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1836 out1file_text=ucase(out1file_text)
1837 if (out1file_text(1:1).eq."Y") then
1840 out1file=fg_rank.gt.0
1845 if (lentmp.gt.0) then
1846 write (*,'(80(1h!))')
1847 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1848 write (*,'(80(1h!))')
1849 write (*,*)"All output files will be on node /tmp directory."
1851 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1852 if (me.eq.king) then
1853 write (*,*) "The master node is ",nodename
1854 else if (fg_rank.eq.0) then
1855 write (*,*) "I am the CG slave node ",nodename
1857 write (*,*) "I am the FG slave node ",nodename
1860 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1861 lenpre = lentmp+lenpre+1
1863 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1864 C Get the names and open the input files
1865 #if defined(WINIFL) || defined(WINPGI)
1866 open(1,file=pref_orig(:ilen(pref_orig))//
1867 & '.inp',status='old',readonly,shared)
1868 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1869 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1870 C Get parameter filenames and open the parameter files.
1871 call getenv_loc('BONDPAR',bondname)
1872 open (ibond,file=bondname,status='old',readonly,shared)
1873 call getenv_loc('THETPAR',thetname)
1874 open (ithep,file=thetname,status='old',readonly,shared)
1875 call getenv_loc('ROTPAR',rotname)
1876 open (irotam,file=rotname,status='old',readonly,shared)
1877 call getenv_loc('TORPAR',torname)
1878 open (itorp,file=torname,status='old',readonly,shared)
1879 call getenv_loc('TORDPAR',tordname)
1880 open (itordp,file=tordname,status='old',readonly,shared)
1881 call getenv_loc('FOURIER',fouriername)
1882 open (ifourier,file=fouriername,status='old',readonly,shared)
1883 call getenv_loc('ELEPAR',elename)
1884 open (ielep,file=elename,status='old',readonly,shared)
1885 call getenv_loc('SIDEPAR',sidename)
1886 open (isidep,file=sidename,status='old',readonly,shared)
1887 #elif (defined CRAY) || (defined AIX)
1888 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1890 c print *,"Processor",myrank," opened file 1"
1891 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1892 c print *,"Processor",myrank," opened file 9"
1893 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1894 C Get parameter filenames and open the parameter files.
1895 call getenv_loc('BONDPAR',bondname)
1896 open (ibond,file=bondname,status='old',action='read')
1897 c print *,"Processor",myrank," opened file IBOND"
1898 call getenv_loc('THETPAR',thetname)
1899 open (ithep,file=thetname,status='old',action='read')
1900 c print *,"Processor",myrank," opened file ITHEP"
1901 call getenv_loc('ROTPAR',rotname)
1902 open (irotam,file=rotname,status='old',action='read')
1903 c print *,"Processor",myrank," opened file IROTAM"
1904 call getenv_loc('TORPAR',torname)
1905 open (itorp,file=torname,status='old',action='read')
1906 c print *,"Processor",myrank," opened file ITORP"
1907 call getenv_loc('TORDPAR',tordname)
1908 open (itordp,file=tordname,status='old',action='read')
1909 c print *,"Processor",myrank," opened file ITORDP"
1910 call getenv_loc('SCCORPAR',sccorname)
1911 open (isccor,file=sccorname,status='old',action='read')
1912 c print *,"Processor",myrank," opened file ISCCOR"
1913 call getenv_loc('FOURIER',fouriername)
1914 open (ifourier,file=fouriername,status='old',action='read')
1915 c print *,"Processor",myrank," opened file IFOURIER"
1916 call getenv_loc('ELEPAR',elename)
1917 open (ielep,file=elename,status='old',action='read')
1918 c print *,"Processor",myrank," opened file IELEP"
1919 call getenv_loc('SIDEPAR',sidename)
1920 open (isidep,file=sidename,status='old',action='read')
1921 c print *,"Processor",myrank," opened file ISIDEP"
1922 c print *,"Processor",myrank," opened parameter files"
1924 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1925 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1926 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1927 C Get parameter filenames and open the parameter files.
1928 call getenv_loc('BONDPAR',bondname)
1929 open (ibond,file=bondname,status='old')
1930 call getenv_loc('THETPAR',thetname)
1931 open (ithep,file=thetname,status='old')
1932 call getenv_loc('ROTPAR',rotname)
1933 open (irotam,file=rotname,status='old')
1934 call getenv_loc('TORPAR',torname)
1935 open (itorp,file=torname,status='old')
1936 call getenv_loc('TORDPAR',tordname)
1937 open (itordp,file=tordname,status='old')
1938 call getenv_loc('SCCORPAR',sccorname)
1939 open (isccor,file=sccorname,status='old')
1940 call getenv_loc('FOURIER',fouriername)
1941 open (ifourier,file=fouriername,status='old')
1942 call getenv_loc('ELEPAR',elename)
1943 open (ielep,file=elename,status='old')
1944 call getenv_loc('SIDEPAR',sidename)
1945 open (isidep,file=sidename,status='old')
1947 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1949 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1950 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1951 C Get parameter filenames and open the parameter files.
1952 call getenv_loc('BONDPAR',bondname)
1953 open (ibond,file=bondname,status='old',readonly)
1954 call getenv_loc('THETPAR',thetname)
1955 open (ithep,file=thetname,status='old',readonly)
1956 call getenv_loc('ROTPAR',rotname)
1957 open (irotam,file=rotname,status='old',readonly)
1958 call getenv_loc('TORPAR',torname)
1959 open (itorp,file=torname,status='old',readonly)
1960 call getenv_loc('TORDPAR',tordname)
1961 open (itordp,file=tordname,status='old',readonly)
1962 call getenv_loc('SCCORPAR',sccorname)
1963 open (isccor,file=sccorname,status='old',readonly)
1965 call getenv_loc('THETPARPDB',thetname_pdb)
1966 print *,"thetname_pdb ",thetname_pdb
1967 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1968 print *,ithep_pdb," opened"
1970 call getenv_loc('FOURIER',fouriername)
1971 open (ifourier,file=fouriername,status='old',readonly)
1972 call getenv_loc('ELEPAR',elename)
1973 open (ielep,file=elename,status='old',readonly)
1974 call getenv_loc('SIDEPAR',sidename)
1975 open (isidep,file=sidename,status='old',readonly)
1977 call getenv_loc('ROTPARPDB',rotname_pdb)
1978 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1983 C 8/9/01 In the newest version SCp interaction constants are read from a file
1984 C Use -DOLDSCP to use hard-coded constants instead.
1986 call getenv_loc('SCPPAR',scpname)
1987 #if defined(WINIFL) || defined(WINPGI)
1988 open (iscpp,file=scpname,status='old',readonly,shared)
1989 #elif (defined CRAY) || (defined AIX)
1990 open (iscpp,file=scpname,status='old',action='read')
1992 open (iscpp,file=scpname,status='old')
1994 open (iscpp,file=scpname,status='old',readonly)
1997 call getenv_loc('PATTERN',patname)
1998 #if defined(WINIFL) || defined(WINPGI)
1999 open (icbase,file=patname,status='old',readonly,shared)
2000 #elif (defined CRAY) || (defined AIX)
2001 open (icbase,file=patname,status='old',action='read')
2003 open (icbase,file=patname,status='old')
2005 open (icbase,file=patname,status='old',readonly)
2008 C Open output file only for CG processes
2009 c print *,"Processor",myrank," fg_rank",fg_rank
2010 if (fg_rank.eq.0) then
2012 if (nodes.eq.1) then
2015 npos = dlog10(dfloat(nodes-1))+1
2017 if (npos.lt.3) npos=3
2018 write (liczba,'(i1)') npos
2019 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2021 write (liczba,form) me
2022 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2023 & liczba(:ilen(liczba))
2024 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2026 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2028 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2029 & liczba(:ilen(liczba))//'.mol2'
2030 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2031 & liczba(:ilen(liczba))//'.stat'
2033 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2034 & //liczba(:ilen(liczba))//'.stat')
2035 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2038 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2039 & liczba(:ilen(liczba))//'.const'
2044 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2045 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2046 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2047 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2048 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2050 & call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
2052 rest2name=prefix(:ilen(prefix))//'.rst'
2054 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2057 #if defined(AIX) || defined(PGI)
2058 if (me.eq.king .or. .not. out1file)
2059 & open(iout,file=outname,status='unknown')
2061 if (fg_rank.gt.0) then
2062 write (liczba,'(i3.3)') myrank/nfgtasks
2063 write (ll,'(bz,i3.3)') fg_rank
2064 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2069 open(igeom,file=intname,status='unknown',position='append')
2070 open(ipdb,file=pdbname,status='unknown')
2071 open(imol2,file=mol2name,status='unknown')
2072 open(istat,file=statname,status='unknown',position='append')
2074 c1out open(iout,file=outname,status='unknown')
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',access='append')
2089 open(ipdb,file=pdbname,status='unknown')
2090 open(imol2,file=mol2name,status='unknown')
2091 open(istat,file=statname,status='unknown',access='append')
2093 c1out open(iout,file=outname,status='unknown')
2096 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2097 csa_seed=prefix(:lenpre)//'.CSA.seed'
2098 csa_history=prefix(:lenpre)//'.CSA.history'
2099 csa_bank=prefix(:lenpre)//'.CSA.bank'
2100 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2101 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2102 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2103 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2104 csa_int=prefix(:lenpre)//'.int'
2105 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2106 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2107 csa_in=prefix(:lenpre)//'.CSA.in'
2108 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2111 write (iout,'(80(1h-))')
2112 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2113 write (iout,'(80(1h-))')
2114 write (iout,*) "Input file : ",
2115 & pref_orig(:ilen(pref_orig))//'.inp'
2116 write (iout,*) "Output file : ",
2117 & outname(:ilen(outname))
2119 write (iout,*) "Sidechain potential file : ",
2120 & sidename(:ilen(sidename))
2122 write (iout,*) "SCp potential file : ",
2123 & scpname(:ilen(scpname))
2125 write (iout,*) "Electrostatic potential file : ",
2126 & elename(:ilen(elename))
2127 write (iout,*) "Cumulant coefficient file : ",
2128 & fouriername(:ilen(fouriername))
2129 write (iout,*) "Torsional parameter file : ",
2130 & torname(:ilen(torname))
2131 write (iout,*) "Double torsional parameter file : ",
2132 & tordname(:ilen(tordname))
2133 write (iout,*) "SCCOR parameter file : ",
2134 & sccorname(:ilen(sccorname))
2135 write (iout,*) "Bond & inertia constant file : ",
2136 & bondname(:ilen(bondname))
2137 write (iout,*) "Bending parameter file : ",
2138 & thetname(:ilen(thetname))
2139 write (iout,*) "Rotamer parameter file : ",
2140 & rotname(:ilen(rotname))
2141 write (iout,*) "Thetpdb parameter file : ",
2142 & thetname_pdb(:ilen(thetname_pdb))
2143 write (iout,*) "Threading database : ",
2144 & patname(:ilen(patname))
2146 &write (iout,*)" DIRTMP : ",
2148 write (iout,'(80(1h-))')
2152 c----------------------------------------------------------------------------
2153 subroutine card_concat(card)
2154 implicit real*8 (a-h,o-z)
2155 include 'DIMENSIONS'
2156 include 'COMMON.IOUNITS'
2158 character*80 karta,ucase
2160 read (inp,'(a)') karta
2163 do while (karta(80:80).eq.'&')
2164 card=card(:ilen(card)+1)//karta(:79)
2165 read (inp,'(a)') karta
2168 card=card(:ilen(card)+1)//karta
2171 c----------------------------------------------------------------------------------
2173 implicit real*8 (a-h,o-z)
2174 include 'DIMENSIONS'
2175 include 'COMMON.CHAIN'
2176 include 'COMMON.IOUNITS'
2178 open(irest2,file=rest2name,status='unknown')
2179 read(irest2,*) totT,EK,potE,totE,t_bath
2181 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2184 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2187 read (irest2,*) iset
2192 c---------------------------------------------------------------------------------
2193 subroutine read_fragments
2194 implicit real*8 (a-h,o-z)
2195 include 'DIMENSIONS'
2199 include 'COMMON.SETUP'
2200 include 'COMMON.CHAIN'
2201 include 'COMMON.IOUNITS'
2203 include 'COMMON.CONTROL'
2204 read(inp,*) nset,nfrag,npair,nfrag_back
2205 if(me.eq.king.or..not.out1file)
2206 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2207 & " nfrag_back",nfrag_back
2209 read(inp,*) mset(iset)
2211 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2213 if(me.eq.king.or..not.out1file)
2214 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2215 & ifrag(2,i,iset), qinfrag(i,iset)
2218 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2220 if(me.eq.king.or..not.out1file)
2221 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2222 & ipair(2,i,iset), qinpair(i,iset)
2225 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2226 & wfrag_back(3,i,iset),
2227 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2228 if(me.eq.king.or..not.out1file)
2229 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2230 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2235 c-------------------------------------------------------------------------------
2236 subroutine read_dist_constr
2237 implicit real*8 (a-h,o-z)
2238 include 'DIMENSIONS'
2242 include 'COMMON.SETUP'
2243 include 'COMMON.CONTROL'
2244 include 'COMMON.CHAIN'
2245 include 'COMMON.IOUNITS'
2246 include 'COMMON.SBRIDGE'
2247 integer ifrag_(2,100),ipair_(2,100)
2248 double precision wfrag_(100),wpair_(100)
2249 character*500 controlcard
2250 c write (iout,*) "Calling read_dist_constr"
2251 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2253 call card_concat(controlcard)
2254 call readi(controlcard,"NFRAG",nfrag_,0)
2255 call readi(controlcard,"NPAIR",npair_,0)
2256 call readi(controlcard,"NDIST",ndist_,0)
2257 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2258 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2259 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2260 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2261 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2262 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2263 c write (iout,*) "IFRAG"
2265 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2267 c write (iout,*) "IPAIR"
2269 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2273 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2274 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2275 & ifrag_(2,i)=nstart_sup+nsup-1
2276 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2278 if (wfrag_(i).gt.0.0d0) then
2279 do j=ifrag_(1,i),ifrag_(2,i)-1
2280 do k=j+1,ifrag_(2,i)
2281 c write (iout,*) "j",j," k",k
2283 if (constr_dist.eq.1) then
2288 forcon(nhpb)=wfrag_(i)
2289 else if (constr_dist.eq.2) then
2290 if (ddjk.le.dist_cut) then
2295 forcon(nhpb)=wfrag_(i)
2302 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2305 if (.not.out1file .or. me.eq.king)
2306 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2307 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2309 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2310 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2317 if (wpair_(i).gt.0.0d0) then
2325 do j=ifrag_(1,ii),ifrag_(2,ii)
2326 do k=ifrag_(1,jj),ifrag_(2,jj)
2330 forcon(nhpb)=wpair_(i)
2331 dhpb(nhpb)=dist(j,k)
2333 if (.not.out1file .or. me.eq.king)
2334 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2335 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2337 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2338 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2345 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2346 if (forcon(nhpb+1).gt.0.0d0) then
2348 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2350 if (.not.out1file .or. me.eq.king)
2351 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2352 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2354 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2355 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2362 c-------------------------------------------------------------------------------
2364 subroutine flush(iu)
2369 subroutine flush(iu)
2374 c------------------------------------------------------------------------------
2375 subroutine copy_to_tmp(source)
2376 include "DIMENSIONS"
2377 include "COMMON.IOUNITS"
2378 character*(*) source
2379 character* 256 tmpfile
2383 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2384 inquire(file=tmpfile,exist=ex)
2386 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2387 & " to temporary directory..."
2388 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2389 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2393 c------------------------------------------------------------------------------
2394 subroutine move_from_tmp(source)
2395 include "DIMENSIONS"
2396 include "COMMON.IOUNITS"
2397 character*(*) source
2400 write (*,*) "Moving ",source(:ilen(source)),
2401 & " from temporary directory to working directory"
2402 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2403 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2406 c------------------------------------------------------------------------------
2407 subroutine random_init(seed)
2409 C Initialize random number generator
2411 implicit real*8 (a-h,o-z)
2412 include 'DIMENSIONS'
2415 logical OKRandom, prng_restart
2417 integer iseed_array(4)
2419 include 'COMMON.IOUNITS'
2420 include 'COMMON.TIME1'
2421 include 'COMMON.THREAD'
2422 include 'COMMON.SBRIDGE'
2423 include 'COMMON.CONTROL'
2424 include 'COMMON.MCM'
2425 include 'COMMON.MAP'
2426 include 'COMMON.HEADER'
2427 include 'COMMON.CSA'
2428 include 'COMMON.CHAIN'
2429 include 'COMMON.MUCA'
2431 include 'COMMON.FFIELD'
2432 include 'COMMON.SETUP'
2433 iseed=-dint(dabs(seed))
2434 if (iseed.eq.0) then
2435 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2436 & 'Random seed undefined. The program will stop.'
2437 write (*,'(/80(1h*)/20x,a/80(1h*))')
2438 & 'Random seed undefined. The program will stop.'
2440 call mpi_finalize(mpi_comm_world,ierr)
2442 stop 'Bad random seed.'
2445 if (fg_rank.eq.0) then
2449 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2450 OKRandom = prng_restart(me,iseed)
2453 tmp=65536.0d0**(4-i)
2454 iseed_array(i) = dint(seed/tmp)
2455 seed=seed-iseed_array(i)*tmp
2458 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2459 & (iseed_array(i),i=1,4)
2460 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2461 & (iseed_array(i),i=1,4)
2462 OKRandom = prng_restart(me,iseed_array)
2465 c r1 = prng_next(me)
2466 r1=ran_number(0.0D0,1.0D0)
2468 & write (iout,*) 'ran_num',r1
2469 if (r1.lt.0.0d0) OKRandom=.false.
2471 if (.not.OKRandom) then
2472 write (iout,*) 'PRNG IS NOT WORKING!!!'
2473 print *,'PRNG IS NOT WORKING!!!'
2476 call mpi_abort(mpi_comm_world,error_msg,ierr)
2479 write (iout,*) 'too many processors for parallel prng'
2480 write (*,*) 'too many processors for parallel prng'
2488 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)