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 if(me.eq.king.or..not.out1file) then
683 write (iout,*) "Parameters of the SS-bond potential:"
684 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
686 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
687 write (iout,*) "EBR",ebr
688 c print *,'indpdb=',indpdb,' pdbref=',pdbref
690 if (indpdb.gt.0 .or. pdbref) then
691 read(inp,'(a)') pdbfile
692 if(me.eq.king.or..not.out1file)
693 & write (iout,'(2a)') 'PDB data will be read from file ',
694 & pdbfile(:ilen(pdbfile))
695 open(ipdbin,file=pdbfile,status='old',err=33)
697 33 write (iout,'(a)') 'Error opening PDB file.'
700 c print *,'Begin reading pdb data'
702 c print *,'Finished reading pdb data'
703 if(me.eq.king.or..not.out1file)
704 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
705 & ' nstart_sup=',nstart_sup
707 itype_pdb(i)=itype(i)
711 nct=nstart_sup+nsup-1
712 call contact(.false.,ncont_ref,icont_ref,co)
715 if(me.eq.king.or..not.out1file)
716 & write(iout,*)'Adding sidechains'
720 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
723 do while (fail.and.nsi.le.maxsi)
724 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
727 if(fail) write(iout,*)'Adding sidechain failed for res ',
728 & i,' after ',nsi,' trials'
733 if (indpdb.eq.0) then
734 C Read sequence if not taken from the pdb file.
736 c print *,'nres=',nres
737 if (iscode.gt.0) then
738 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
740 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
742 C Convert sequence to numeric code
744 itype(i)=rescode(i,sequence(i),iscode)
746 C Assign initial virtual bond lengths
752 vbld(i+nres)=dsc(iabs(itype(i)))
753 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
754 c write (iout,*) "i",i," itype",itype(i),
755 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
759 c print '(20i4)',(itype(i),i=1,nres)
762 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
764 if (itype(i).eq.ntyp1) then
768 else if (iabs(itype(i+1)).ne.20) then
770 else if (iabs(itype(i)).ne.20) then
777 if(me.eq.king.or..not.out1file)then
778 write (iout,*) "ITEL"
780 write (iout,*) i,itype(i),itel(i)
782 print *,'Call Read_Bridge.'
785 C 8/13/98 Set limits to generating the dihedral angles
790 read (inp,*) ndih_constr
791 if (ndih_constr.gt.0) then
793 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
794 if(me.eq.king.or..not.out1file)then
796 & 'There are',ndih_constr,' constraints on phi angles.'
798 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
802 phi0(i)=deg2rad*phi0(i)
803 drange(i)=deg2rad*drange(i)
805 if(me.eq.king.or..not.out1file)
806 & write (iout,*) 'FTORS',ftors
809 phibound(1,ii) = phi0(i)-drange(i)
810 phibound(2,ii) = phi0(i)+drange(i)
817 write (iout,'(a)') 'Boundaries in phi angle sampling:'
819 write (iout,'(a3,i5,2f10.1)')
820 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
826 cd print *,'NNT=',NNT,' NCT=',NCT
827 if (itype(1).eq.ntyp1) nnt=2
828 if (itype(nres).eq.ntyp1) nct=nct-1
830 if(me.eq.king.or..not.out1file)
831 & write (iout,'(a,i3)') 'nsup=',nsup
833 if (nsup.le.(nct-nnt+1)) then
834 do i=0,nct-nnt+1-nsup
835 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
841 & 'Error - sequences to be superposed do not match.'
844 do i=0,nsup-(nct-nnt+1)
845 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
847 nstart_sup=nstart_sup+i
853 & 'Error - sequences to be superposed do not match.'
856 if (nsup.eq.0) nsup=nct-nnt
857 if (nstart_sup.eq.0) nstart_sup=nnt
858 if (nstart_seq.eq.0) nstart_seq=nnt
859 if(me.eq.king.or..not.out1file)
860 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
861 & ' nstart_seq=',nstart_seq
863 c--- Zscore rms -------
864 if (nz_start.eq.0) nz_start=nnt
865 if (nz_end.eq.0 .and. nsup.gt.0) then
867 else if (nz_end.eq.0) then
870 if(me.eq.king.or..not.out1file)then
871 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
872 write (iout,*) 'IZ_SC=',iz_sc
874 c----------------------
877 if (.not.pdbref) then
878 call read_angles(inp,*38)
880 38 write (iout,'(a)') 'Error reading reference structure.'
882 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
883 stop 'Error reading reference structure'
887 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
897 call contact(.true.,ncont_ref,icont_ref,co)
899 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
901 if (constr_dist.gt.0) call read_dist_constr
902 write (iout,*) "After read_dist_constr nhpb",nhpb
904 if(me.eq.king.or..not.out1file)
905 & write (iout,*) 'Contact order:',co
907 if(me.eq.king.or..not.out1file)
908 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
911 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
913 if(me.eq.king.or..not.out1file)
914 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
915 & icont_ref(1,i),' ',
916 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
920 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
921 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
922 & modecalc.ne.10) then
923 C If input structure hasn't been supplied from the PDB file read or generate
925 if (iranconf.eq.0 .and. .not. extconf) then
926 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
927 & write (iout,'(a)') 'Initial geometry will be read in.'
929 read(inp,'(8f10.5)',end=36,err=36)
930 & ((c(l,k),l=1,3),k=1,nres),
931 & ((c(l,k+nres),l=1,3),k=nnt,nct)
932 write (iout,*) "Exit READ_CART"
933 write (iout,'(8f10.5)')
934 & ((c(l,k),l=1,3),k=1,nres),
935 & ((c(l,k+nres),l=1,3),k=nnt,nct)
936 call int_from_cart1(.true.)
937 write (iout,*) "Finish INT_TO_CART"
940 dc(j,i)=c(j,i+1)-c(j,i)
941 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
945 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
947 dc(j,i+nres)=c(j,i+nres)-c(j,i)
948 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
954 call read_angles(inp,*36)
957 36 write (iout,'(a)') 'Error reading angle file.'
959 call mpi_finalize( MPI_COMM_WORLD,IERR )
961 stop 'Error reading angle file.'
963 else if (extconf) then
964 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
965 & write (iout,'(a)') 'Extended chain initial geometry.'
967 theta(i)=90d0*deg2rad
973 alph(i)=110d0*deg2rad
976 omeg(i)=-120d0*deg2rad
977 if (itype(i).le.0) omeg(i)=-omeg(i)
980 if(me.eq.king.or..not.out1file)
981 & write (iout,'(a)') 'Random-generated initial geometry.'
985 if (me.eq.king .or. fg_rank.eq.0 .and. (
986 & modecalc.eq.12 .or. modecalc.eq.14) ) then
990 call gen_rand_conf(itmp,*30)
992 30 write (iout,*) 'Failed to generate random conformation',
994 write (*,*) 'Processor:',me,
995 & ' Failed to generate random conformation',
1005 write (iout,'(a,i3,a)') 'Processor:',me,
1006 & ' error in generating random conformation.'
1007 write (*,'(a,i3,a)') 'Processor:',me,
1008 & ' error in generating random conformation.'
1011 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1017 call gen_rand_conf(itmp,*30)
1019 30 write (iout,*) 'Failed to generate random conformation',
1020 & ', itrial=',itrial
1021 write (*,*) 'Failed to generate random conformation',
1022 & ', itrial=',itrial
1024 write (iout,'(a,i3,a)') 'Processor:',me,
1025 & ' error in generating random conformation.'
1026 write (*,'(a,i3,a)') 'Processor:',me,
1027 & ' error in generating random conformation.'
1032 elseif (modecalc.eq.4) then
1033 read (inp,'(a)') intinname
1034 open (intin,file=intinname,status='old',err=333)
1035 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1036 & write (iout,'(a)') 'intinname',intinname
1037 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1039 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1041 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1043 stop 'Error opening angle file.'
1047 C Generate distance constraints, if the PDB structure is to be regularized.
1048 if (nthread.gt.0) then
1049 call read_threadbase
1052 if (me.eq.king .or. .not. out1file)
1054 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1055 write (iout,'(/a,i3,a)')
1056 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1057 write (iout,'(20i4)') (iss(i),i=1,ns)
1058 write (iout,'(/a/)') 'Pre-formed links are:'
1064 if (me.eq.king.or..not.out1file)
1065 & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1066 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1071 if (i2ndstr.gt.0) call secstrp2dihc
1072 c call geom_to_var(nvar,x)
1073 c call etotal(energia(0))
1074 c call enerprint(energia(0))
1075 c call briefout(0,etot)
1077 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1078 cd write (iout,'(a)') 'Variable list:'
1079 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1081 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1082 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1083 & 'Processor',myrank,': end reading molecular data.'
1087 c--------------------------------------------------------------------------
1088 logical function seq_comp(itypea,itypeb,length)
1090 integer length,itypea(length),itypeb(length)
1093 if (itypea(i).ne.itypeb(i)) then
1101 c-----------------------------------------------------------------------------
1102 subroutine read_bridge
1103 C Read information about disulfide bridges.
1104 implicit real*8 (a-h,o-z)
1105 include 'DIMENSIONS'
1109 include 'COMMON.IOUNITS'
1110 include 'COMMON.GEO'
1111 include 'COMMON.VAR'
1112 include 'COMMON.INTERACT'
1113 include 'COMMON.LOCAL'
1114 include 'COMMON.NAMES'
1115 include 'COMMON.CHAIN'
1116 include 'COMMON.FFIELD'
1117 include 'COMMON.SBRIDGE'
1118 include 'COMMON.HEADER'
1119 include 'COMMON.CONTROL'
1120 include 'COMMON.DBASE'
1121 include 'COMMON.THREAD'
1122 include 'COMMON.TIME1'
1123 include 'COMMON.SETUP'
1124 C Read bridging residues.
1125 read (inp,*) ns,(iss(i),i=1,ns)
1127 if(me.eq.king.or..not.out1file)
1128 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1129 C Check whether the specified bridging residues are cystines.
1131 if (itype(iss(i)).ne.1) then
1132 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1133 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
1134 & ' can form a disulfide bridge?!!!'
1135 write (*,'(2a,i3,a)')
1136 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
1137 & ' can form a disulfide bridge?!!!'
1139 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1144 C Read preformed bridges.
1146 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1147 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1150 C Check if the residues involved in bridges are in the specified list of
1151 C bridging residues.
1154 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1155 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1156 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1157 & ' contains residues present in other pairs.'
1158 write (*,'(a,i3,a)') 'Disulfide pair',i,
1159 & ' contains residues present in other pairs.'
1161 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1167 if (ihpb(i).eq.iss(j)) goto 10
1169 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1172 if (jhpb(i).eq.iss(j)) goto 20
1174 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1180 ihpb(i)=ihpb(i)+nres
1181 jhpb(i)=jhpb(i)+nres
1187 c----------------------------------------------------------------------------
1188 subroutine read_x(kanal,*)
1189 implicit real*8 (a-h,o-z)
1190 include 'DIMENSIONS'
1191 include 'COMMON.GEO'
1192 include 'COMMON.VAR'
1193 include 'COMMON.CHAIN'
1194 include 'COMMON.IOUNITS'
1195 include 'COMMON.CONTROL'
1196 include 'COMMON.LOCAL'
1197 include 'COMMON.INTERACT'
1198 c Read coordinates from input
1200 read(kanal,'(8f10.5)',end=10,err=10)
1201 & ((c(l,k),l=1,3),k=1,nres),
1202 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1205 c(j,2*nres)=c(j,nres)
1207 call int_from_cart1(.false.)
1210 dc(j,i)=c(j,i+1)-c(j,i)
1211 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1215 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1217 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1218 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1226 c----------------------------------------------------------------------------
1227 subroutine read_threadbase
1228 implicit real*8 (a-h,o-z)
1229 include 'DIMENSIONS'
1230 include 'COMMON.IOUNITS'
1231 include 'COMMON.GEO'
1232 include 'COMMON.VAR'
1233 include 'COMMON.INTERACT'
1234 include 'COMMON.LOCAL'
1235 include 'COMMON.NAMES'
1236 include 'COMMON.CHAIN'
1237 include 'COMMON.FFIELD'
1238 include 'COMMON.SBRIDGE'
1239 include 'COMMON.HEADER'
1240 include 'COMMON.CONTROL'
1241 include 'COMMON.DBASE'
1242 include 'COMMON.THREAD'
1243 include 'COMMON.TIME1'
1244 C Read pattern database for threading.
1245 read (icbase,*) nseq
1247 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1248 & nres_base(2,i),nres_base(3,i)
1249 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1251 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1252 c & nres_base(2,i),nres_base(3,i)
1253 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1257 if (weidis.eq.0.0D0) weidis=0.1D0
1266 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1267 write (iout,'(a,i5)') 'nexcl: ',nexcl
1268 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1271 c------------------------------------------------------------------------------
1272 subroutine setup_var
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 Set up variable list.
1295 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1297 ialph(i,1)=nvar+nside
1301 if (indphi.gt.0) then
1303 else if (indback.gt.0) then
1308 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1311 c----------------------------------------------------------------------------
1312 subroutine gen_dist_constr
1313 C Generate CA distance constraints.
1314 implicit real*8 (a-h,o-z)
1315 include 'DIMENSIONS'
1316 include 'COMMON.IOUNITS'
1317 include 'COMMON.GEO'
1318 include 'COMMON.VAR'
1319 include 'COMMON.INTERACT'
1320 include 'COMMON.LOCAL'
1321 include 'COMMON.NAMES'
1322 include 'COMMON.CHAIN'
1323 include 'COMMON.FFIELD'
1324 include 'COMMON.SBRIDGE'
1325 include 'COMMON.HEADER'
1326 include 'COMMON.CONTROL'
1327 include 'COMMON.DBASE'
1328 include 'COMMON.THREAD'
1329 include 'COMMON.TIME1'
1330 dimension itype_pdb(maxres)
1331 common /pizda/ itype_pdb
1333 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1334 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1335 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1337 do i=nstart_sup,nstart_sup+nsup-1
1338 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1339 cd & ' seq_pdb', restyp(itype_pdb(i))
1340 do j=i+2,nstart_sup+nsup-1
1342 ihpb(nhpb)=i+nstart_seq-nstart_sup
1343 jhpb(nhpb)=j+nstart_seq-nstart_sup
1345 dhpb(nhpb)=dist(i,j)
1348 cd write (iout,'(a)') 'Distance constraints:'
1353 cd if (ii.gt.nres) then
1358 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1359 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1360 cd & dhpb(i),forcon(i)
1364 c----------------------------------------------------------------------------
1366 implicit real*8 (a-h,o-z)
1367 include 'DIMENSIONS'
1368 include 'COMMON.MAP'
1369 include 'COMMON.IOUNITS'
1370 character*3 angid(4) /'THE','PHI','ALP','OME'/
1371 character*80 mapcard,ucase
1373 read (inp,'(a)') mapcard
1374 mapcard=ucase(mapcard)
1375 if (index(mapcard,'PHI').gt.0) then
1377 else if (index(mapcard,'THE').gt.0) then
1379 else if (index(mapcard,'ALP').gt.0) then
1381 else if (index(mapcard,'OME').gt.0) then
1384 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1385 stop 'Error - illegal variable spec in MAP card.'
1387 call readi (mapcard,'RES1',res1(imap),0)
1388 call readi (mapcard,'RES2',res2(imap),0)
1389 if (res1(imap).eq.0) then
1390 res1(imap)=res2(imap)
1391 else if (res2(imap).eq.0) then
1392 res2(imap)=res1(imap)
1394 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1396 & 'Error - illegal definition of variable group in MAP.'
1397 stop 'Error - illegal definition of variable group in MAP.'
1399 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1400 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1401 call readi(mapcard,'NSTEP',nstep(imap),0)
1402 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1404 & 'Illegal boundary and/or step size specification in MAP.'
1405 stop 'Illegal boundary and/or step size specification in MAP.'
1410 c----------------------------------------------------------------------------
1412 implicit real*8 (a-h,o-z)
1413 include 'DIMENSIONS'
1414 include 'COMMON.IOUNITS'
1415 include 'COMMON.GEO'
1416 include 'COMMON.CSA'
1417 include 'COMMON.BANK'
1418 include 'COMMON.CONTROL'
1420 character*620 mcmcard
1421 call card_concat(mcmcard)
1423 call readi(mcmcard,'NCONF',nconf,50)
1424 call readi(mcmcard,'NADD',nadd,0)
1425 call readi(mcmcard,'JSTART',jstart,1)
1426 call readi(mcmcard,'JEND',jend,1)
1427 call readi(mcmcard,'NSTMAX',nstmax,500000)
1428 call readi(mcmcard,'N0',n0,1)
1429 call readi(mcmcard,'N1',n1,6)
1430 call readi(mcmcard,'N2',n2,4)
1431 call readi(mcmcard,'N3',n3,0)
1432 call readi(mcmcard,'N4',n4,0)
1433 call readi(mcmcard,'N5',n5,0)
1434 call readi(mcmcard,'N6',n6,10)
1435 call readi(mcmcard,'N7',n7,0)
1436 call readi(mcmcard,'N8',n8,0)
1437 call readi(mcmcard,'N9',n9,0)
1438 call readi(mcmcard,'N14',n14,0)
1439 call readi(mcmcard,'N15',n15,0)
1440 call readi(mcmcard,'N16',n16,0)
1441 call readi(mcmcard,'N17',n17,0)
1442 call readi(mcmcard,'N18',n18,0)
1444 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1446 call readi(mcmcard,'NDIFF',ndiff,2)
1447 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1448 call readi(mcmcard,'IS1',is1,1)
1449 call readi(mcmcard,'IS2',is2,8)
1450 call readi(mcmcard,'NRAN0',nran0,4)
1451 call readi(mcmcard,'NRAN1',nran1,2)
1452 call readi(mcmcard,'IRR',irr,1)
1453 call readi(mcmcard,'NSEED',nseed,20)
1454 call readi(mcmcard,'NTOTAL',ntotal,10000)
1455 call reada(mcmcard,'CUT1',cut1,2.0d0)
1456 call reada(mcmcard,'CUT2',cut2,5.0d0)
1457 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1458 call readi(mcmcard,'ICMAX',icmax,3)
1459 call readi(mcmcard,'IRESTART',irestart,0)
1460 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1463 call reada(mcmcard,'DELE',dele,20.0d0)
1464 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1465 call readi(mcmcard,'IREF',iref,0)
1466 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1467 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1468 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1469 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1470 write (iout,*) "NCONF_IN",nconf_in
1473 c----------------------------------------------------------------------------
1474 cfmc subroutine mcmfread
1475 cfmc implicit real*8 (a-h,o-z)
1476 cfmc include 'DIMENSIONS'
1477 cfmc include 'COMMON.MCMF'
1478 cfmc include 'COMMON.IOUNITS'
1479 cfmc include 'COMMON.GEO'
1480 cfmc character*80 ucase
1481 cfmc character*620 mcmcard
1482 cfmc call card_concat(mcmcard)
1484 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1485 cfmc write(iout,*)'MAXRANT=',maxrant
1486 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1487 cfmc write(iout,*)'MAXFAM=',maxfam
1488 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1489 cfmc write(iout,*)'NNET1=',nnet1
1490 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1491 cfmc write(iout,*)'NNET2=',nnet2
1492 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1493 cfmc write(iout,*)'NNET3=',nnet3
1494 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1495 cfmc write(iout,*)'ILASTT=',ilastt
1496 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1497 cfmc write(iout,*)'MAXSTR=',maxstr
1498 cfmc maxstr_f=maxstr/maxfam
1499 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1500 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1501 cfmc write(iout,*)'NMCMF=',nmcmf
1502 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1503 cfmc write(iout,*)'IFOCUS=',ifocus
1504 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1505 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1506 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1507 cfmc write(iout,*)'INTPRT=',intprt
1508 cfmc call readi(mcmcard,'IPRT',iprt,100)
1509 cfmc write(iout,*)'IPRT=',iprt
1510 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1511 cfmc write(iout,*)'IMAXTR=',imaxtr
1512 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1513 cfmc write(iout,*)'MAXEVEN=',maxeven
1514 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1515 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1516 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1517 cfmc write(iout,*)'INIMIN=',inimin
1518 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1519 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1520 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1521 cfmc write(iout,*)'NTHREAD=',nthread
1522 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1523 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1524 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1525 cfmc write(iout,*)'MAXPERT=',maxpert
1526 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1527 cfmc write(iout,*)'IRMSD=',irmsd
1528 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1529 cfmc write(iout,*)'DENEMIN=',denemin
1530 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1531 cfmc write(iout,*)'RCUT1S=',rcut1s
1532 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1533 cfmc write(iout,*)'RCUT1E=',rcut1e
1534 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1535 cfmc write(iout,*)'RCUT2S=',rcut2s
1536 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1537 cfmc write(iout,*)'RCUT2E=',rcut2e
1538 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1539 cfmc write(iout,*)'DPERT1=',d_pert1
1540 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1541 cfmc write(iout,*)'DPERT1A=',d_pert1a
1542 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1543 cfmc write(iout,*)'DPERT2=',d_pert2
1544 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1545 cfmc write(iout,*)'DPERT2A=',d_pert2a
1546 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1547 cfmc write(iout,*)'DPERT2B=',d_pert2b
1548 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1549 cfmc write(iout,*)'DPERT2C=',d_pert2c
1550 cfmc d_pert1=deg2rad*d_pert1
1551 cfmc d_pert1a=deg2rad*d_pert1a
1552 cfmc d_pert2=deg2rad*d_pert2
1553 cfmc d_pert2a=deg2rad*d_pert2a
1554 cfmc d_pert2b=deg2rad*d_pert2b
1555 cfmc d_pert2c=deg2rad*d_pert2c
1556 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1557 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1558 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1559 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1560 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1561 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1562 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1563 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1564 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1565 cfmc write(iout,*)'RCUTINI=',rcutini
1566 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1567 cfmc write(iout,*)'GRAT=',grat
1568 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1569 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1573 c----------------------------------------------------------------------------
1575 implicit real*8 (a-h,o-z)
1576 include 'DIMENSIONS'
1577 include 'COMMON.MCM'
1578 include 'COMMON.MCE'
1579 include 'COMMON.IOUNITS'
1581 character*320 mcmcard
1582 call card_concat(mcmcard)
1583 call readi(mcmcard,'MAXACC',maxacc,100)
1584 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1585 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1586 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1587 call readi(mcmcard,'MAXREPM',maxrepm,200)
1588 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1589 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1590 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1591 call reada(mcmcard,'E_UP',e_up,5.0D0)
1592 call reada(mcmcard,'DELTE',delte,0.1D0)
1593 call readi(mcmcard,'NSWEEP',nsweep,5)
1594 call readi(mcmcard,'NSTEPH',nsteph,0)
1595 call readi(mcmcard,'NSTEPC',nstepc,0)
1596 call reada(mcmcard,'TMIN',tmin,298.0D0)
1597 call reada(mcmcard,'TMAX',tmax,298.0D0)
1598 call readi(mcmcard,'NWINDOW',nwindow,0)
1599 call readi(mcmcard,'PRINT_MC',print_mc,0)
1600 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1601 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1602 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1603 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1604 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1605 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1606 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1607 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1608 if (nwindow.gt.0) then
1609 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1611 winlen(i)=winend(i)-winstart(i)+1
1614 if (tmax.lt.tmin) tmax=tmin
1615 if (tmax.eq.tmin) then
1619 if (nstepc.gt.0 .and. nsteph.gt.0) then
1620 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1621 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1623 C Probabilities of different move types
1624 sumpro_type(0)=0.0D0
1625 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1626 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1627 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1628 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1629 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1630 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1631 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1633 print *,'i',i,' sumprotype',sumpro_type(i)
1634 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1635 print *,'i',i,' sumprotype',sumpro_type(i)
1639 c----------------------------------------------------------------------------
1640 subroutine read_minim
1641 implicit real*8 (a-h,o-z)
1642 include 'DIMENSIONS'
1643 include 'COMMON.MINIM'
1644 include 'COMMON.IOUNITS'
1646 character*320 minimcard
1647 call card_concat(minimcard)
1648 call readi(minimcard,'MAXMIN',maxmin,2000)
1649 call readi(minimcard,'MAXFUN',maxfun,5000)
1650 call readi(minimcard,'MINMIN',minmin,maxmin)
1651 call readi(minimcard,'MINFUN',minfun,maxmin)
1652 call reada(minimcard,'TOLF',tolf,1.0D-2)
1653 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1654 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1655 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1656 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1657 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1658 & 'Options in energy minimization:'
1659 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1660 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1661 & 'MinMin:',MinMin,' MinFun:',MinFun,
1662 & ' TolF:',TolF,' RTolF:',RTolF
1665 c----------------------------------------------------------------------------
1666 subroutine read_angles(kanal,*)
1667 implicit real*8 (a-h,o-z)
1668 include 'DIMENSIONS'
1669 include 'COMMON.GEO'
1670 include 'COMMON.VAR'
1671 include 'COMMON.CHAIN'
1672 include 'COMMON.IOUNITS'
1673 include 'COMMON.CONTROL'
1674 c Read angles from input
1676 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1677 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1678 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1679 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1682 c 9/7/01 avoid 180 deg valence angle
1683 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1685 theta(i)=deg2rad*theta(i)
1686 phi(i)=deg2rad*phi(i)
1687 alph(i)=deg2rad*alph(i)
1688 omeg(i)=deg2rad*omeg(i)
1693 c----------------------------------------------------------------------------
1694 subroutine reada(rekord,lancuch,wartosc,default)
1696 character*(*) rekord,lancuch
1697 double precision wartosc,default
1700 iread=index(rekord,lancuch)
1701 if (iread.eq.0) then
1705 iread=iread+ilen(lancuch)+1
1706 read (rekord(iread:),*,err=10,end=10) wartosc
1711 c----------------------------------------------------------------------------
1712 subroutine readi(rekord,lancuch,wartosc,default)
1714 character*(*) rekord,lancuch
1715 integer wartosc,default
1718 iread=index(rekord,lancuch)
1719 if (iread.eq.0) then
1723 iread=iread+ilen(lancuch)+1
1724 read (rekord(iread:),*,err=10,end=10) wartosc
1729 c----------------------------------------------------------------------------
1730 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1733 integer tablica(dim),default
1734 character*(*) rekord,lancuch
1741 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1742 if (iread.eq.0) return
1743 iread=iread+ilen(lancuch)+1
1744 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1747 c----------------------------------------------------------------------------
1748 subroutine multreada(rekord,lancuch,tablica,dim,default)
1751 double precision tablica(dim),default
1752 character*(*) rekord,lancuch
1759 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1760 if (iread.eq.0) return
1761 iread=iread+ilen(lancuch)+1
1762 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1765 c----------------------------------------------------------------------------
1766 subroutine openunits
1767 implicit real*8 (a-h,o-z)
1768 include 'DIMENSIONS'
1771 character*16 form,nodename
1774 include 'COMMON.SETUP'
1775 include 'COMMON.IOUNITS'
1777 include 'COMMON.CONTROL'
1778 integer lenpre,lenpot,ilen,lentmp
1780 character*3 out1file_text,ucase
1783 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1784 call getenv_loc("PREFIX",prefix)
1786 call getenv_loc("POT",pot)
1787 call getenv_loc("DIRTMP",tmpdir)
1788 call getenv_loc("CURDIR",curdir)
1789 call getenv_loc("OUT1FILE",out1file_text)
1790 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1791 out1file_text=ucase(out1file_text)
1792 if (out1file_text(1:1).eq."Y") then
1795 out1file=fg_rank.gt.0
1800 if (lentmp.gt.0) then
1801 write (*,'(80(1h!))')
1802 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1803 write (*,'(80(1h!))')
1804 write (*,*)"All output files will be on node /tmp directory."
1806 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1807 if (me.eq.king) then
1808 write (*,*) "The master node is ",nodename
1809 else if (fg_rank.eq.0) then
1810 write (*,*) "I am the CG slave node ",nodename
1812 write (*,*) "I am the FG slave node ",nodename
1815 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1816 lenpre = lentmp+lenpre+1
1818 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1819 C Get the names and open the input files
1820 #if defined(WINIFL) || defined(WINPGI)
1821 open(1,file=pref_orig(:ilen(pref_orig))//
1822 & '.inp',status='old',readonly,shared)
1823 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1824 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1825 C Get parameter filenames and open the parameter files.
1826 call getenv_loc('BONDPAR',bondname)
1827 open (ibond,file=bondname,status='old',readonly,shared)
1828 call getenv_loc('THETPAR',thetname)
1829 open (ithep,file=thetname,status='old',readonly,shared)
1830 call getenv_loc('ROTPAR',rotname)
1831 open (irotam,file=rotname,status='old',readonly,shared)
1832 call getenv_loc('TORPAR',torname)
1833 open (itorp,file=torname,status='old',readonly,shared)
1834 call getenv_loc('TORDPAR',tordname)
1835 open (itordp,file=tordname,status='old',readonly,shared)
1836 call getenv_loc('FOURIER',fouriername)
1837 open (ifourier,file=fouriername,status='old',readonly,shared)
1838 call getenv_loc('ELEPAR',elename)
1839 open (ielep,file=elename,status='old',readonly,shared)
1840 call getenv_loc('SIDEPAR',sidename)
1841 open (isidep,file=sidename,status='old',readonly,shared)
1842 #elif (defined CRAY) || (defined AIX)
1843 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1845 c print *,"Processor",myrank," opened file 1"
1846 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1847 c print *,"Processor",myrank," opened file 9"
1848 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1849 C Get parameter filenames and open the parameter files.
1850 call getenv_loc('BONDPAR',bondname)
1851 open (ibond,file=bondname,status='old',action='read')
1852 c print *,"Processor",myrank," opened file IBOND"
1853 call getenv_loc('THETPAR',thetname)
1854 open (ithep,file=thetname,status='old',action='read')
1855 c print *,"Processor",myrank," opened file ITHEP"
1856 call getenv_loc('ROTPAR',rotname)
1857 open (irotam,file=rotname,status='old',action='read')
1858 c print *,"Processor",myrank," opened file IROTAM"
1859 call getenv_loc('TORPAR',torname)
1860 open (itorp,file=torname,status='old',action='read')
1861 c print *,"Processor",myrank," opened file ITORP"
1862 call getenv_loc('TORDPAR',tordname)
1863 open (itordp,file=tordname,status='old',action='read')
1864 c print *,"Processor",myrank," opened file ITORDP"
1865 call getenv_loc('SCCORPAR',sccorname)
1866 open (isccor,file=sccorname,status='old',action='read')
1867 c print *,"Processor",myrank," opened file ISCCOR"
1868 call getenv_loc('FOURIER',fouriername)
1869 open (ifourier,file=fouriername,status='old',action='read')
1870 c print *,"Processor",myrank," opened file IFOURIER"
1871 call getenv_loc('ELEPAR',elename)
1872 open (ielep,file=elename,status='old',action='read')
1873 c print *,"Processor",myrank," opened file IELEP"
1874 call getenv_loc('SIDEPAR',sidename)
1875 open (isidep,file=sidename,status='old',action='read')
1876 c print *,"Processor",myrank," opened file ISIDEP"
1877 c print *,"Processor",myrank," opened parameter files"
1879 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1880 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1881 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1882 C Get parameter filenames and open the parameter files.
1883 call getenv_loc('BONDPAR',bondname)
1884 open (ibond,file=bondname,status='old')
1885 call getenv_loc('THETPAR',thetname)
1886 open (ithep,file=thetname,status='old')
1887 call getenv_loc('ROTPAR',rotname)
1888 open (irotam,file=rotname,status='old')
1889 call getenv_loc('TORPAR',torname)
1890 open (itorp,file=torname,status='old')
1891 call getenv_loc('TORDPAR',tordname)
1892 open (itordp,file=tordname,status='old')
1893 call getenv_loc('SCCORPAR',sccorname)
1894 open (isccor,file=sccorname,status='old')
1895 call getenv_loc('FOURIER',fouriername)
1896 open (ifourier,file=fouriername,status='old')
1897 call getenv_loc('ELEPAR',elename)
1898 open (ielep,file=elename,status='old')
1899 call getenv_loc('SIDEPAR',sidename)
1900 open (isidep,file=sidename,status='old')
1902 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1904 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1905 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1906 C Get parameter filenames and open the parameter files.
1907 call getenv_loc('BONDPAR',bondname)
1908 open (ibond,file=bondname,status='old',readonly)
1909 call getenv_loc('THETPAR',thetname)
1910 open (ithep,file=thetname,status='old',readonly)
1911 call getenv_loc('ROTPAR',rotname)
1912 open (irotam,file=rotname,status='old',readonly)
1913 call getenv_loc('TORPAR',torname)
1914 open (itorp,file=torname,status='old',readonly)
1915 call getenv_loc('TORDPAR',tordname)
1916 open (itordp,file=tordname,status='old',readonly)
1917 call getenv_loc('SCCORPAR',sccorname)
1918 open (isccor,file=sccorname,status='old',readonly)
1919 call getenv_loc('FOURIER',fouriername)
1920 open (ifourier,file=fouriername,status='old',readonly)
1921 call getenv_loc('ELEPAR',elename)
1922 open (ielep,file=elename,status='old',readonly)
1923 call getenv_loc('SIDEPAR',sidename)
1924 open (isidep,file=sidename,status='old',readonly)
1928 C 8/9/01 In the newest version SCp interaction constants are read from a file
1929 C Use -DOLDSCP to use hard-coded constants instead.
1931 call getenv_loc('SCPPAR',scpname)
1932 #if defined(WINIFL) || defined(WINPGI)
1933 open (iscpp,file=scpname,status='old',readonly,shared)
1934 #elif (defined CRAY) || (defined AIX)
1935 open (iscpp,file=scpname,status='old',action='read')
1937 open (iscpp,file=scpname,status='old')
1939 open (iscpp,file=scpname,status='old',readonly)
1942 call getenv_loc('PATTERN',patname)
1943 #if defined(WINIFL) || defined(WINPGI)
1944 open (icbase,file=patname,status='old',readonly,shared)
1945 #elif (defined CRAY) || (defined AIX)
1946 open (icbase,file=patname,status='old',action='read')
1948 open (icbase,file=patname,status='old')
1950 open (icbase,file=patname,status='old',readonly)
1953 C Open output file only for CG processes
1954 c print *,"Processor",myrank," fg_rank",fg_rank
1955 if (fg_rank.eq.0) then
1957 if (nodes.eq.1) then
1960 npos = dlog10(dfloat(nodes-1))+1
1962 if (npos.lt.3) npos=3
1963 write (liczba,'(i1)') npos
1964 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1966 write (liczba,form) me
1967 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1968 & liczba(:ilen(liczba))
1969 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1971 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1973 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1974 & liczba(:ilen(liczba))//'.mol2'
1975 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1976 & liczba(:ilen(liczba))//'.stat'
1978 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1979 & //liczba(:ilen(liczba))//'.stat')
1980 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1983 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1984 & liczba(:ilen(liczba))//'.const'
1989 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1990 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1991 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1992 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1993 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1995 & call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
1997 rest2name=prefix(:ilen(prefix))//'.rst'
1999 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2002 #if defined(AIX) || defined(PGI)
2003 if (me.eq.king .or. .not. out1file)
2004 & open(iout,file=outname,status='unknown')
2006 if (fg_rank.gt.0) then
2007 write (liczba,'(i3.3)') myrank/nfgtasks
2008 write (ll,'(bz,i3.3)') fg_rank
2009 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2014 open(igeom,file=intname,status='unknown',position='append')
2015 open(ipdb,file=pdbname,status='unknown')
2016 open(imol2,file=mol2name,status='unknown')
2017 open(istat,file=statname,status='unknown',position='append')
2019 c1out open(iout,file=outname,status='unknown')
2022 if (me.eq.king .or. .not.out1file)
2023 & open(iout,file=outname,status='unknown')
2025 if (fg_rank.gt.0) then
2026 write (liczba,'(i3.3)') myrank/nfgtasks
2027 write (ll,'(bz,i3.3)') fg_rank
2028 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2033 open(igeom,file=intname,status='unknown',access='append')
2034 open(ipdb,file=pdbname,status='unknown')
2035 open(imol2,file=mol2name,status='unknown')
2036 open(istat,file=statname,status='unknown',access='append')
2038 c1out open(iout,file=outname,status='unknown')
2041 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2042 csa_seed=prefix(:lenpre)//'.CSA.seed'
2043 csa_history=prefix(:lenpre)//'.CSA.history'
2044 csa_bank=prefix(:lenpre)//'.CSA.bank'
2045 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2046 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2047 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2048 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2049 csa_int=prefix(:lenpre)//'.int'
2050 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2051 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2052 csa_in=prefix(:lenpre)//'.CSA.in'
2053 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2056 write (iout,'(80(1h-))')
2057 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2058 write (iout,'(80(1h-))')
2059 write (iout,*) "Input file : ",
2060 & pref_orig(:ilen(pref_orig))//'.inp'
2061 write (iout,*) "Output file : ",
2062 & outname(:ilen(outname))
2064 write (iout,*) "Sidechain potential file : ",
2065 & sidename(:ilen(sidename))
2067 write (iout,*) "SCp potential file : ",
2068 & scpname(:ilen(scpname))
2070 write (iout,*) "Electrostatic potential file : ",
2071 & elename(:ilen(elename))
2072 write (iout,*) "Cumulant coefficient file : ",
2073 & fouriername(:ilen(fouriername))
2074 write (iout,*) "Torsional parameter file : ",
2075 & torname(:ilen(torname))
2076 write (iout,*) "Double torsional parameter file : ",
2077 & tordname(:ilen(tordname))
2078 write (iout,*) "SCCOR parameter file : ",
2079 & sccorname(:ilen(sccorname))
2080 write (iout,*) "Bond & inertia constant file : ",
2081 & bondname(:ilen(bondname))
2082 write (iout,*) "Bending parameter file : ",
2083 & thetname(:ilen(thetname))
2084 write (iout,*) "Rotamer parameter file : ",
2085 & rotname(:ilen(rotname))
2086 write (iout,*) "Threading database : ",
2087 & patname(:ilen(patname))
2089 &write (iout,*)" DIRTMP : ",
2091 write (iout,'(80(1h-))')
2095 c----------------------------------------------------------------------------
2096 subroutine card_concat(card)
2097 implicit real*8 (a-h,o-z)
2098 include 'DIMENSIONS'
2099 include 'COMMON.IOUNITS'
2101 character*80 karta,ucase
2103 read (inp,'(a)') karta
2106 do while (karta(80:80).eq.'&')
2107 card=card(:ilen(card)+1)//karta(:79)
2108 read (inp,'(a)') karta
2111 card=card(:ilen(card)+1)//karta
2114 c----------------------------------------------------------------------------------
2116 implicit real*8 (a-h,o-z)
2117 include 'DIMENSIONS'
2118 include 'COMMON.CHAIN'
2119 include 'COMMON.IOUNITS'
2121 open(irest2,file=rest2name,status='unknown')
2122 read(irest2,*) totT,EK,potE,totE,t_bath
2124 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2127 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2130 read (irest2,*) iset
2135 c---------------------------------------------------------------------------------
2136 subroutine read_fragments
2137 implicit real*8 (a-h,o-z)
2138 include 'DIMENSIONS'
2142 include 'COMMON.SETUP'
2143 include 'COMMON.CHAIN'
2144 include 'COMMON.IOUNITS'
2146 include 'COMMON.CONTROL'
2147 read(inp,*) nset,nfrag,npair,nfrag_back
2148 if(me.eq.king.or..not.out1file)
2149 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2150 & " nfrag_back",nfrag_back
2152 read(inp,*) mset(iset)
2154 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2156 if(me.eq.king.or..not.out1file)
2157 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2158 & ifrag(2,i,iset), qinfrag(i,iset)
2161 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2163 if(me.eq.king.or..not.out1file)
2164 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2165 & ipair(2,i,iset), qinpair(i,iset)
2168 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2169 & wfrag_back(3,i,iset),
2170 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2171 if(me.eq.king.or..not.out1file)
2172 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2173 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2178 c-------------------------------------------------------------------------------
2179 subroutine read_dist_constr
2180 implicit real*8 (a-h,o-z)
2181 include 'DIMENSIONS'
2185 include 'COMMON.SETUP'
2186 include 'COMMON.CONTROL'
2187 include 'COMMON.CHAIN'
2188 include 'COMMON.IOUNITS'
2189 include 'COMMON.SBRIDGE'
2190 integer ifrag_(2,100),ipair_(2,100)
2191 double precision wfrag_(100),wpair_(100)
2192 character*500 controlcard
2193 c write (iout,*) "Calling read_dist_constr"
2194 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2196 call card_concat(controlcard)
2197 call readi(controlcard,"NFRAG",nfrag_,0)
2198 call readi(controlcard,"NPAIR",npair_,0)
2199 call readi(controlcard,"NDIST",ndist_,0)
2200 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2201 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2202 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2203 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2204 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2205 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2206 c write (iout,*) "IFRAG"
2208 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2210 c write (iout,*) "IPAIR"
2212 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2216 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2217 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2218 & ifrag_(2,i)=nstart_sup+nsup-1
2219 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2221 if (wfrag_(i).gt.0.0d0) then
2222 do j=ifrag_(1,i),ifrag_(2,i)-1
2223 do k=j+1,ifrag_(2,i)
2224 c write (iout,*) "j",j," k",k
2226 if (constr_dist.eq.1) then
2231 forcon(nhpb)=wfrag_(i)
2232 else if (constr_dist.eq.2) then
2233 if (ddjk.le.dist_cut) then
2238 forcon(nhpb)=wfrag_(i)
2245 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2248 if (.not.out1file .or. me.eq.king)
2249 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2250 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2252 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2253 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2260 if (wpair_(i).gt.0.0d0) then
2268 do j=ifrag_(1,ii),ifrag_(2,ii)
2269 do k=ifrag_(1,jj),ifrag_(2,jj)
2273 forcon(nhpb)=wpair_(i)
2274 dhpb(nhpb)=dist(j,k)
2276 if (.not.out1file .or. me.eq.king)
2277 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2278 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2280 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2281 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2288 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2289 if (forcon(nhpb+1).gt.0.0d0) then
2291 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2293 if (.not.out1file .or. me.eq.king)
2294 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2295 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2297 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2298 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2305 c-------------------------------------------------------------------------------
2307 subroutine flush(iu)
2312 subroutine flush(iu)
2317 c------------------------------------------------------------------------------
2318 subroutine copy_to_tmp(source)
2319 include "DIMENSIONS"
2320 include "COMMON.IOUNITS"
2321 character*(*) source
2322 character* 256 tmpfile
2326 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2327 inquire(file=tmpfile,exist=ex)
2329 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2330 & " to temporary directory..."
2331 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2332 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2336 c------------------------------------------------------------------------------
2337 subroutine move_from_tmp(source)
2338 include "DIMENSIONS"
2339 include "COMMON.IOUNITS"
2340 character*(*) source
2343 write (*,*) "Moving ",source(:ilen(source)),
2344 & " from temporary directory to working directory"
2345 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2346 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2349 c------------------------------------------------------------------------------
2350 subroutine random_init(seed)
2352 C Initialize random number generator
2354 implicit real*8 (a-h,o-z)
2355 include 'DIMENSIONS'
2358 logical OKRandom, prng_restart
2360 integer iseed_array(4)
2362 include 'COMMON.IOUNITS'
2363 include 'COMMON.TIME1'
2364 include 'COMMON.THREAD'
2365 include 'COMMON.SBRIDGE'
2366 include 'COMMON.CONTROL'
2367 include 'COMMON.MCM'
2368 include 'COMMON.MAP'
2369 include 'COMMON.HEADER'
2370 include 'COMMON.CSA'
2371 include 'COMMON.CHAIN'
2372 include 'COMMON.MUCA'
2374 include 'COMMON.FFIELD'
2375 include 'COMMON.SETUP'
2376 iseed=-dint(dabs(seed))
2377 if (iseed.eq.0) then
2378 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2379 & 'Random seed undefined. The program will stop.'
2380 write (*,'(/80(1h*)/20x,a/80(1h*))')
2381 & 'Random seed undefined. The program will stop.'
2383 call mpi_finalize(mpi_comm_world,ierr)
2385 stop 'Bad random seed.'
2388 if (fg_rank.eq.0) then
2392 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2393 OKRandom = prng_restart(me,iseed)
2396 tmp=65536.0d0**(4-i)
2397 iseed_array(i) = dint(seed/tmp)
2398 seed=seed-iseed_array(i)*tmp
2401 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2402 & (iseed_array(i),i=1,4)
2403 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2404 & (iseed_array(i),i=1,4)
2405 OKRandom = prng_restart(me,iseed_array)
2408 c r1 = prng_next(me)
2409 r1=ran_number(0.0D0,1.0D0)
2411 & write (iout,*) 'ran_num',r1
2412 if (r1.lt.0.0d0) OKRandom=.false.
2414 if (.not.OKRandom) then
2415 write (iout,*) 'PRNG IS NOT WORKING!!!'
2416 print *,'PRNG IS NOT WORKING!!!'
2419 call mpi_abort(mpi_comm_world,error_msg,ierr)
2422 write (iout,*) 'too many processors for parallel prng'
2423 write (*,*) 'too many processors for parallel prng'
2431 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)